{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Learning Tree-augmented Naive Bayes (TAN) Structure from Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we show an example for learning the structure of a Bayesian Network using the TAN algorithm.  We will first build a model to generate some data and then attempt to learn the model's graph structure back from the generated data.\n",
    "\n",
    "For comparison of Naive Bayes and TAN classifier, refer to the blog post [Classification with TAN and Pgmpy](https://loudly-soft.blogspot.com/2020/08/classification-with-tree-augmented.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First, create a Naive Bayes graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAjYElEQVR4nO3dWW8bZ4Lu8YdVXGVR+2rJsS3ZjuN4TU4HSeckGDjoJJNz05gGBg3M3MzlfItzNXfzBQ7QfTd9c3qAIL0hpweZdM900t3psZ04TuzYlmNrMSVKJkVKIllkVZ0LjRXJNiWK4lKs+v+AIF64lKVX78Oqp+qtkOu6rgAACAij3RsAAEArEXwAgEAh+AAAgULwAQACheADAAQKwQcACBSCDwAQKAQfACBQCD4AQKAQfACAQCH4AACBQvABAAKF4AMABArBBwAIFIIPABAoBB8AIFDC7d4AeIvjuLJsR64rhUJS1DRkGKF2bxawL4xj7Ibgg7IbluYyBaVWC0qvWbJdd+vvzFBIw91RjfUmNNmfUF9XtI1bClTHOEatQq67bXQgUB6uFnR9blWLuaJMI6SuWFiJiClz2ydj23FVKNvaKFVkO65Ge+I6N9mr8d5EG7cc+A7jGPtF8AVQsWzr6mxGt1NrSibC6o6FFQrtfRjIdV2tlSrKFco6NZbUpSP9ikfMFmwx8DTGMepF8AVMdsPSRzeXVCo7GkxGZdQwUTzJcV2t5C3FIoYunx7hsBFajnGMgyD4AiS7YenDrxYVNUNKxiMHfr18sSzLdvXOmVEmDbQM4xgHxeUMAVEs2/ro5lLDJgtJSsYjipohfXRzScWy3ZDXBHbDOEYjEHwBcXU2o1LZadhk8VgyHlGp7OjqbKahrws8C+MYjUDwBcDD1YJup9Y0mGzOYZzBZFTfpPJKrRab8vqAxDhG4xB8AXB9blXJRLiuEwBqYYRC6klE9MVctimvD0iMYzQOF7D7XHbD0mKuqLHe+J6PvfLJ71UqbX7aDYWkcCSqZE+vjp54XvFE167P7Y6FlVotKrthcYIAGu6g47i3f0BHTzyvaDS263MZx8HAHp/PzWUKMo1QTdc3PdY/OKTRiSMyzbAeLac1c+urPZ8TCoVkGiHNZQoH2Vzgmeoex4ePSJKWF1O6f+fWns9hHAcDwedzqdWCumL727EfGZ/U8ZMv6NiJ5yVJhY31mp7XFQsrtcqEgcarexyfekHPTZ2UJK3nczU9j3Hsfxzq9DHHcZVeszRwaH+HbJYezmk1u6LVRyuSpIHhkZqel4iYSq9ZchyXBYHRMPWO483nOlpf2wy8ru5kTc9hHPsfwedjlu3Idt0daxbWIrOyvPVrwwjpUHdPTc8zjZBsd3NV/LjBElBojHrH8a0vr239uqevX8dPvlDT8xjH/kfw+Vi9a/I8f/aiBoZHlF/N6qtrn2nm1g319PXveYLLQd8XeJZ6x1P/4JAc19XqoxWt5bIqbKwrEq19r5Fx7F90fD520LO+u3t6ZRimXFcqFjZa9r7AdvWOp5HxSZ258LKGxw7LcVx9e+dmS94X3scen49FTUNmKCTb2d9hoscd33o+r0qlIsMwaupHbMeVGQopavJ5Co1T7zh+7MjxaS0vLmg9n1f20bL6BoZ2fTzj2P/4zvqYYWzefLOwz/UHMyvLSs3NqrC+pmRvn06fv7Tn9U+SVCjbGu6OckIAGqrecfxYLJ7Q0OhhSdL8/Xt7Pp5x7H/cncHnvpxf1RdzWY307H3h70Et5Yo6P9mnsxO9NT3etm1VKhWZpqlwmIMPqM7L4xidh9nG5yb7E7r6ICPXdfd18e9+ua6rxXRa93P3VFg4pHK5rFKppHK5LMuytv57/PtKpSLXdWXbti5evKjXX3+9aduGztfKcWw7rib7uTO7nxF8PtfXFdVoT1xrpUrDV7Tfbq1UUU/E1a/f/7+KRqM6fPiwBgcHFY1GZRiGTNOUaZqKRqNKJpMyDEOhUEizs7M6fvx407YL/tDKcTzaE2e5Mp+j4wuAc5O9yhcqcpp0VNtxXeUKZf2vV8/pH//xHzU6OipJmp2dVT6fVzweV1dXl+LxuCKRiEzTVCgU0tramgYGBjQ+Pt6U7YK/tGocn5/sa8rrwzsIvgAY703o5Fi3VvJWU15/JW/p1FhSY71xvfjii3r99dfV3d2tF154QZVKRV988YVmZma0sbHzkohMJqOXX365qYeu4B+tHMfwN4IvIC4d6VcsYihfLDf0dfPFsmIRQ5eO9EvaXOT3tdde0+HDh5XL5XTs2DGdP39e8Xhct27d0s2bN5XNZlUsFhWPx3Xs2LGGbg/8rVXjGP5G8AVEPGLq8ukRWbbbsEkjXyzLsl1dPj2ieOS7pZ3C4bB+8IMfKBaLaXV1VZFIRIcPH9aFCxc0NDSkubk5ffLJJxofH5dhMARRu1aOY/gXs06A9HVF9c6ZUbmulM6V6u5KHNdVOleS60rvnBl95okAXV1deu+997S+vq5icfPeaIZhaGhoSC+88IImJyd16NAh/du//Ztu3bqlUql0oH8bgqPR4/j+/Qd6/blDnNASIFzHF0DFsq2rsxl9k8qrJxFRdyxcU8/muq7WShXlCmWdGkvq0pH+PT8h37t3T7/61a80MTGxda1eKpXS6dOn9cYbb2htbU337t3T/Py8xsbGNDU1pZ6e2hbFRrA1ahx/8sHPtLKU0g9/+EOdPXtWXV21rUmLzkXwBVhqtagv5rJazBVlGiF1xcJKRMwdy0LZjqtC2dZGqSLbcTXaE9f5yb59nQBw5coVffrppzpyZPOmoHNzc/rxj3+sgYGBrcdYlqX79+/r22+/VTKZ1NTUlIaHhznxBXs66Dj+wx/+oH/913/VyMiIBgYGdPHiRQLQ5wg+KLthaS5TUGq1oPSaJXvbkDBDm8tFjfUmNNmfqOtwkOM4+uijj3T37l0lEgkNDQ3pvffeq/rYhYUF3b17V47jaGpqSpOTkzJNuhfsrt5xPDMzow8++ECrq6s6efKk8vm8QqEQAehjBB92cJzN+5C57ubq9FHTaMiahZZl6Re/+IXu37+vH//4xzp8+PCez1leXtbMzIwymYyOHTumY8eOKRbbe81QYD/jeHl5WT//+c8VjUaVSqV05swZua6rpaUlAtCnCD60TD6f1/Xr1/Xaa6/t6xAmPSCaybIs/fSnP9XExITu3r0rwzC2VhOqVCoEoA8RfOgY9IBoln/5l39RLBZTOBzWjRs3dPjwYQ0NfXf7oicD8Hvf+x7jroNxOQM6RjQa1cmTJ/XWW29pcnJSX3/9tT7++GPdv39ftl3fLWsASRodHdXGxoZM09SJEyf04MEDFQqFrb8Ph8MaHBxUqVTa8efoTOzxoaPRA6IRvvzyy61FFSQpnU5v9X2maSqTyahQKOjy5cs6ceIEe3sdjuCDL9AD4iDm5+f1y1/+csdJV3fv3pUkJRIJ9fb26u23395xCQ46F4c64Qvd3d06d+6cLl++rEOHDulPf/qT/vjHP2ppaUl8tsNeenp6nhon4+PjmpmZ0cDAgP7mb/6G0PMR7scHX3ncA05PT2thYUFff/21bty4wfWA2NWhQ4dkmqZs295xaPMf/uEftLi4KMuyFI2ypJlfcKgTvkcPiFq8//772tjYUD6f33Fo88GDB5qZmdEbb7zBByefIPgQGPSA2M0f//hH/f73v9crr7yi73//+zv28K5evSrDMHThwoU2biEaheBD4HA9IJ4lnU4rl8tpamrqqbFQqVT0H//xHzp58qQmJyfbtIVoFIIPgcW6oNiPXC6nTz/9VK+//rq6u7vbvTk4AIIPED0gakPf5w8EH7ANPSD2cuXKFZmmSd/XwQg+4BnoAVENfV/nI/iAXdAD4lno+zobwQfUiB4Q29H3dS6CD9gnekA8Rt/XmQg+oE70gKDv60wEH3BA9IDBRt/XeQg+oIHoAYOJvq+zEHxAE9ADBg99X+cg+IAmogcMDvq+zkHwAS1ADxgM9H2dgeADWowe0N/o+7yP4APahB7Qv+j7vI3gA9qMHtB/6Pu8jeADPIIe0F/o+7yL4AM8iB7QH+j7vIngAzyMHrDz0fd5D8EHdAB6wM5F3+c9BB/QQegBOxN9n7cQfECHogfsLPR93kHwAR2OHrBz0Pd5A8EH+AQ9oPfR93kDwQf4DD2gt9H3tR/BB/gYPaA30fe1F8EHBAA9oPfQ97UPwQcECD2gd9D3tQ/BBwQQPaA30Pe1B8EHBBw9YHvR97UewQdAEj1gO9H3tRbBB2AHesDWo+9rLYIPwDPRA7YWfV/rEHwA9kQP2Br0fa1B8AGoGT1g89H3NR/BB2Df6AGbh76v+Qg+AHWjB2wO+r7mIvgANAQ9YGPR9zUPwQegoegBG4e+rzkIPgBNQQ94cPR9zUHwAWgqesCDoe9rPIIPQMvQA9aHvq+xCD4ALUcPuH/0fY1D8AFoG3rA2tH3NQ7BB6Dt6AFrQ9/XGAQfAE+hB9wdfd/BEXwAPIkesDr6voMh+AB4Gj3g0+j7DobgA9AR6AF3ou+rH8EHoOPQA26i76sPwQegY9ED0vfVg+AD0PGC3APS9+0fwQfAN4LaA9L37Q/BB8CXgtYD0vfVjuAD4GtB6gHp+2pD8AEIhCD0gPR9tSH4AASK33tA+r69EXwAAsuvPSB93+4IPgCB58cekL6vOoIPAP6bn3pA+r7qCD4AeIJfekD6vmcj+ABgF53eA9L3PY3gA4AadHIPSN+3E8EHAPvQiT0gfd9OBB8A1KHTekD6vu8QfABwQNt7wKNHj+r48eOe7AHp+zYRfADQIJ3QA169elWGYQS67yP4AKDBtveA3d3dmp6e9kwPSN9H8AFA03i1Bwx630fwAUALeK0HDHLfR/ABQAt5qQcMat9H8AFAG3ihBwxq30fwAUAbtbsHDGLfR/ABgEe0qwcMWt9H8AGAx7SjBwxS30fwAYBHtbIHDFLfR/ABgMe1qgcMSt9H8AFAB2l2DxiEvq+jgs9xXFm2I9eVQiEpahoyjPYvAQQArdbMHrDWvq9T52TPB192w9JcpqDUakHpNUv2ts01QyENd0c11pvQZH9CfV3RNm4pALReM3rA3fo+P8zJng2+h6sFXZ9b1WKuKNMIqSsWViJiytz2acJ2XBXKtjZKFdmOq9GeuM5N9mq8N9HGLQeA1mt0D/hk3+enOdlzwVcs27o6m9Ht1JqSibC6Y+GaPrm4rqu1UkW5QlmnxpK6dKRf8Yg/j08DwG5WVlZ09+7dA/eAs7OzWs5kVewe99Wc7Kngy25Y+ujmkkplR4PJqIw6dtUd19VK3lIsYujy6RHP7moDQLOtr69rZmam7h7Qr3OyZ4Ivu2Hpw68WFTVDSsYjB369fLEsy3b1zplRT3yhAaBd6ukB/TwneyL4imVbv77+UKGQGvIFfixfLMt1pffOjXtmFxsA2qXWHtDvc7Ingu/TmWV9m97QcE/j16RL50o6Ntyl16aGGv7aANCpdusB/T4ntz34Hq4W9Nsbixrri9d1/HgvjuvqYbagd14c11hvvObnWZal27dva3x8XAMDAw3fLgDwgid7wIHJaX38zYrn5uRGanvw/b8bKa1blYbuTj8pXyzrUDSst18c2/OxlmXp1q1b+vOf/6yVlRX96Ec/0vT0dNO2DQC8wLIszc/P66tcWJbTmF6vmv3Myc0Qbsu7/rfshqXFXLHm1M9lM1p4cE/51awcx1YkGlf/0JCOTj8vwzCqPq87FlZqtajshlW1VN0eeJZlaXh4WOVyua5/FwB0mmg0qv7RCT1KLdQ8J1/55PcqlYqSpIuvvK7EoUM1Pa+WObmZ2hp8c5mCTCNU0zUhy4sPdefr63Jdqau7W93JXpVKRS3Oz+rIsRO7Bl8oFJJphDSXKTz1RX5W4LXi/lcA4DX7mZNz2UdboSdJ6cUFPTd1sqb32W1OboW2Bl9qtaCu2N6bYNu27n3ztVxXGh4d1/QLZ7e+MYWNdRk1rEzQFQsrtVrQ2YleSQQeADyp1jlZktKph5KkQ8mk1vN5Lace1hx80tNzciu1Lfgcx1V6zdLAob3TPr+aUaVSkSRNHJ3a8Wkk0VXbrnUiYiq9ZqlYLOn27W8IPADYZj9zsuM4epRelCQdnT6lb778XKVSUbnsI/X01XYy4OM52XHcli9s3bbgs2xHtuvuWOetmsq2ri2WqG/NN9MIaXZuXv/7/f8ju1zS4OCgIpGI7t+/X/U5S0tL+vzzz7WyslLXewJAp3BCpmx3qKY5ObOSVqVSUSQaVU/fgPoHh5VefKh06mHNwWcaIdnu5t0d4kZrr+lrW/Dt51zScOS7s4tKhULNBeqT+vv7dfLV13T75g0ZhqHBwcFdF3C1bVsTExM6evRoXe8HAJ2i7IZ0536hpscupxYkSf2Dm6u/DAyPKr34UI/Sizp+6oVdz7l4UjuuK2hb8O3n8pBkb7/C4bAqlYrm78/s6PhKxYIi0VhNX+hEV0J//T/f1dtv/ZWuX7+uzz//XKFQSCMjI88MwPX1dQ0NDWl0dLT2jQWADlQs29L9uT0fVymXlX20LElaejivpYfz3/1dpaLMSlqDw7XPmU24VHBPbQu+qGnIDIVkO3sf7jRNU8dPvaA7X19XevGh1tfz6k72qmyVtJpZ0cvf/6s9g892XJmh0OaNEg8d0quvvqpz587VFIAA4He1zskrSyk5jivTNNXb/91hzY31dRULG1pOLdQUfNvn5FZrW/AZxuYNC9csW901nEU0NDquaCym+fvfai2X1fLGuiLRuEbGJ2s6q7NQtjXcHd1Roh4iAAFAUu1zcnpx82zO0cOTOnri+a0/z2Uf6cbVvyj7aFlly1IkuvtJMs+ak1ulrSu3fDm/qi/mshrpaf6yNUu5os5P9u166uz6+vqOACyXy3r33XdZuQVAIHhtTm6W1u9jbjPZn5DtuGp29rquK9txNdm/+xmhj/cA//7v/17nz5+XYRg1XcgJAJ2sUqloZmZG3964oo1iyTNzcrOwVucuNjY2FIlEFIk0b9sAoF0KhYLu3bunBw8eaHh4WFNTU/psoeTZOblR2rpyiySdm+zVb28s6lAs3LSVwHOFcl23wOjq6mr49gBAu2WzWd29e1fpdFpHjhzRm2++uTXfnTMKnp2TG6XtwTfem9DJse6m3ftpJW/p1Fiybbe/AAAvcF1XqVRKMzMzKhaLOn78uC5cuKBweGcMBGFObnvwSdKlI/16mC0qXyw3/G6/sYihS0f6G/aaANBJKpWKHjx4oHv37ikWi2l6elpjY2O7nr/g9zm57R3fY9kNSx9+taio2Zj7QOWLZVm2q3fOjLZl9W8AaKdn9Xf9/bUHTnbD0m+uL8ixyxrt7znw9nhpTvZM8EmbX+iPbi6pVHY0mIzWdXzZcV2t5C3FIoYunx5p+xcYAFrpyf7u+PHjdZ+v8NEf/qSvsiENDo35ak72VPBJm8vmXJ3N6JtUXj2JiLpj4ZouKXBdV2ulinKFsk6NJXXpSL/iES5CB+B/z+rvnnvuuaf6u/2YmZnR/Py8/scrr+ra/Kqv5mTPBd9jqdWivpjLajFXlGmE1BULKxExdyylYzuuCmVbG6WKbMfVaE9c5yf7OJEFQCDU09/VIpPJ6M9//rPeeOONrb1FP83Jng2+x7IbluYyBaVWC0qvWbK3ba4Z2lxiZ6w3ocn+hCd2oQGg2Q7a3+2mXC7rd7/7nc6ePauxsaevs/PDnOz54NvOcTbv3eS6myt6R02jLeu8AUA7NLK/q+azzz5TV1eXXnzxxT0f26lzsicuZ6iVYYRafsNCAGinWq+/a4TH7/Hyyy/X9PhOnZM7KvgAICia1d9Vk8lkdPv2bb3xxhv7upFsJyL4AMBDnuzvXnrppYb1d9WUy2X913/9ly5cuBCIpRo7quMDAL9qRX9XzX56PT9gjw8A2qSV/V01++31/IDgA4AWa3V/V02Qer3tCD4AaJF29HfVBK3X246ODwCarJ39XTVB6/W2Y48PAJrAC/1dNUHs9bZr/3cAAHzEK/1dNUHt9bYj+ACgAbzU31UT5F5vOzo+ADgAL/Z31QS519uOPT4A2Ccv93fVBL3X28673yUA8Biv93fV0OvtRPABwB46ob+rhl7vaXR8AFBFJ/V31dDrPY09PgDYphP7u2ro9Z6t876TANAEndrfVUOvVx3BByDQOrm/q4Zeb3d0fAACyQ/9XTX0ertjjw9AYPipv6uGXm9v/vluA0AVfuvvqqHXqw3BB8C3/NjfVUOvVzs6PgC+4+f+rhp6vdqxxwfAF4LQ31VDr7c//h8RAHwtKP1dNfR6+0fwAehIQervqqHXqw8dH4COEsT+rhp6vfqwxwfA84Lc31VDr1e/4I4aAJ4X9P6uGnq9gyH4AHgO/V119HoHR8cHwDPo7/ZGr3dw7PEBaCv6u9rR6zUGIwtAW9Df7Q+9XuMQfABaiv5u/+j1GouOD0BL0N/Vj16vsdjjA9A09HcHR6/XeIw+AA1Hf9cY9HrNQfABaBj6u8ah12seOj4AB0Z/13j0es3DHh+AutDfNQ+9XnMxQgHsC/1dc9HrNR/BB6Am9HfNR6/XGnR8AHZFf9c69HqtwR4fgKfQ37UevV7rMIoBbKG/aw96vdYi+ADQ37URvV7r0fEBAUZ/1370eq3HHh8QMPR33kGv1x6MdCAg6O+8hV6vfQg+wOfo77yHXq+96PgAn6K/8y56vfZijw/wEfo776PXaz9+GgAfoL/rDPR63kDwAR2M/q5z0Ot5Bx0f0IHo7zoPvZ53sMcHdAj6u85Fr+ct/MQAHkd/19no9byH4AM8iv6u89HreRMdH+Ax9Hf+Qa/nTezxAR5Af+c/9HrexU8V0Eb0d/5Er+dtBB/QBvR3/kWv5310fEAL0d/5H72e97HHBzQZ/V1w0Ot1Bn7ygCahvwsWer3OQfABDUZ/Fzz0ep2Fjg9oEPq74KLX6yzs8QEHQH8Her3Ow08nUIft/V08HtfU1BT9XQDR63Umgg/YB/o7PEav17no+IAa0N/hSfR6nYs9PqAK+jtUQ6/X2fgJBp5Af4fd0Ot1PoIP+G/0d9gLvZ4/0PEh8OjvUCt6PX9gjw+BRH+HJz3eB6h2SJtezz/4KUeg0N+hmm+++UY3b97UW2+9pe7u7h1/R6/nLxzqRCA82d9NTU3R32GHf//3f9fnn3+uZDKpt99+W0eOHJG02ev97ne/09mzZzU2NtbmrUQjEHzwNfo71OpnP/uZIpGIHMfR8vKyXnnlFb300ku6cuUKvZ7PEHzwnWf1d8899xz9Haoql8v6yU9+osOHDysUCsm2bc3PzyscDuvYsWP6wQ9+wCFOH2EmgG/Q36FeuVxO0ncntpimqYGBAX322WdyHEfz8/Nbhz7R+Qg+dDyuv8NB5XI5bT/4ValUdOfOHV28eFGxWEwffPDB1qFP0zTbuKVoBIIPHevJ/u7NN9+kv0Nd0um0IpHI1u/v3bungYGBrQ9QExMT+uyzz7SwsPDMsz7RWQg+dBSuv0MzpFKprQ9NqVRKlmVpenp66+9N09TExIRu374ty7L0ox/9iM6vgzFboKWKxaLi8fi+n0d/h2ZxXVdLS0saGhrS2tqaFhYW9OKLL24F2+OzPC3L0ksvvaRLly4Reh2O4EPL3LlzR5988on+9m//tubwo79Ds62vr6tSqch1Xd25c0fHjx9XLBbbEXjPP/+8Ll26xNjzCYIPLZFOp/Xb3/5W5XJZd+/e3fOaqEwmo5mZGfo7NN3jMzof93q9vb1aWloi8HyM4MMOjuPKsh25rhQKSVHTkGEc7HDi+vq6fv3rX6u3t1eRSER/+ctfdPr06afOjqO/Q6PsZxxns1ktLy/LMAwlk0ktLCwQeD7HjAJlNyzNZQpKrRaUXrNkbzut2wyFNNwd1VhvQpP9CfV1Rff12uVyWR9++KFs29bAwICkzb2/2dlZHTt2TBL9HRqj3nF8+/ZtPXjwQGfOnNH09DSBFwCs3BJgD1cLuj63qsVcUaYRUlcsrETElLntk7HtuCqUbW2UKrIdV6M9cZ2b7NV4b2LP13ddVx9//LFu3bqliYmJrT/P5XKKxWJ677339O2337J+Jg7koOP4n//5n3X06FG99dZbjL+AIPgCqFi2dXU2o9upNSUTYXXHwjXtXbmuq7VSRblCWafGkrp0pF/xSPWLea9du6b//M//1HPPPbfj9dfW1nTt2jVdvHhR58+fZ/1M1KVR4/jC4R51xfd3JAOdjeALmOyGpY9uLqlUdjSYjMqo43Ci47payVuKRQxdPj3yzMOf9+/f1y9+8QtNTEwoHA7LdV1lMhmlUimVy2VFo1GdP39e7777biP+WQiYVo1j+BPBFyDZDUsffrWoqBlSMh7Z+wl7yBfLsmxX75wZ3TFpPHr0SD//+c/V19enaDSqdDqtxcVFRSIRjY2Nqb+/X67ramFhQX/3d3+nnp6eA28LgqNV4xj+xVWYAVEs2/ro5lLDJgtJSsYjipohfXRzScWyLWnzurvf/OY3Mk1T6XRa165d09ramqanp3XmzBkNDAwoFArJMAyZpqmvvvqqIduCYGjVOIa/EXwBcXU2o1LZadhk8VgyHlGp7OjqbEa2bev999/XzZs3NTs7K0k6e/asTpw48cy1DYeHh3X9+nUVi8WGbhP8qxXjGP7H5QwB8HC1oNupNY317X+psFoMJqP6JpXX/NfX9Omnn2pkZEQDAwMyTVMrKys7Vr13XVemaW7t8eVyuZouaAdaNY6PD3ZrrLc57wFvIPgC4PrcqpKJcF0nANTCCIXUk4go0XNG//RPr8pxHFUqFdm2/dT/y+WyLMva+q9UKnGbF9SkVeP4i7msxnrHmvIe8AaCz+eyG5YWc8U9P8Fe+eT3KpU2DzkahqFIJKru3l6NTx5Vsrdvz/fpjoWVWi2qYkTV180JAmisWsextHMsb3f+e6/pUHdy1+c+HsfZDYsTXXyM4PO5uUxBphGqeRWU/sEhhSNR5VezWlla1KP0ok6cOa+hkd0/AYdCIZlGSHOZAhMGGm6/41jaHMvxxHfXh26/3141jONgIPh8LrVaUFes9m/zyPikBoZHNleq//q6lhdTunfrK/UPDu95SLIrFlZqtaCzE70H3Wxgh/2OY+m7sbxfjGP/I/h8zHFcpdcsDRza/yfXUCikyaPTWl5MqVKpKL+aUd/A0K7PSURMpdcsOY574IWtgcfqHcdLD+eUyz7a+v2xk6dreh7j2P8IPh+zbEe26+5Ys3A/otvumVcpl/d8vGmEZLubq+LHDU5YQWPUO44zK8s7fl9r8DGO/Y/g87GDrsljbbu+LlxDP9Ko9wW2q3c8PX/2Yl2HOg/6vvA+LmD3sYOc9e26rubu35UkhcNhJXtrX7Weuwmhkdo1nhjH/sUen49FTUNmKCTbqf0w0dLDOT1aXlJ+NatiYUOhkHT8+TM1XWtnO67MUEhRk89TaJx6xrH0dMc3PD6x5+UMEuM4CAg+HzOMzZtvrlm2ums8Iy6zsrx1Hd/gyKgOHzmm7p7azm4rlG0Nd0c5IQANVc84lp7u+Hr6BmoKPsax/xF8PjfWm9AXc9k9J4yXvv/mgd9ro1TRieGn1+QEDqrWcSwdfCwzjv2PfXmfm+xPyHZcNfvuU67rynZcTfbvfWd2YL8Yx2gkgs/n+rqiGu2Ja61Uaer7rJUqGu2Js9oFmoJxjEYi+ALg3GSv8oWKnCZ9WnZcV7lCWecn+5ry+oDEOEbjEHwBMN6b0Mmxbq3kraa8/kre0qmxJLdyQVMxjtEoBF9AXDrSr1jEUL649wos+5EvlhWLGLp0pPbr/IB6MY7RCARfQMQjpi6fHpFluw2bNPLFsizb1eXTI4pHWNoJzcc4RiMQfAHS1xXVO2dG5bpSOlequytxXFfpXEmuK71zZpQTAdBSjGMcVMht9vnB8Jxi2dbV2Yy+SeXVk4ioOxau6T5nrutqrVRRrlDWqbGkLh3p5xMy2oZxjHoRfAGWWi3qi7msFnNFmUZIXbGwEhFzx7JQtuOqULa1UarIdlyN9sR1frKPEwDgGYxj7BfBB2U3LM1lCkqtFpRes2RvGxJmaHO5qLHehCb7ExwOgmcxjlErgg87OM7mfchcd3N1+qhpsGYhOg7jGLsh+AAAgcJZnQCAQCH4AACBQvABAAKF4AMABArBBwAIFIIPABAoBB8AIFAIPgBAoBB8AIBAIfgAAIFC8AEAAoXgAwAECsEHAAgUgg8AECgEHwAgUP4/otz6IArFaOkAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import networkx as nx\n",
    "import matplotlib.pyplot as plt\n",
    "from pgmpy.models import DiscreteBayesianNetwork\n",
    "\n",
    "# class variable is A and feature variables are B, C, D, E and R\n",
    "model = DiscreteBayesianNetwork([(\"A\", \"R\"), (\"A\", \"B\"), (\"A\", \"C\"), (\"A\", \"D\"), (\"A\", \"E\")])\n",
    "nx.draw_circular(\n",
    "    model, with_labels=True, arrowsize=30, node_size=800, alpha=0.3, font_weight=\"bold\"\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Second, add interaction between the features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAv7klEQVR4nO3d6XMbd37n8U934yQJgOAJXpJ46KJE6rK8tjVSYnlHdryZJFXZSk0mebCpfbRPt2qr9uE+StXuf7D7fCaVrCeZ8kxmop2JJ75GsSWLOiyZkniI4oWLBAhCOBro7n3AkCJFgriBPj6vKpdtkSAh6cf+Av3uQ9A0TQMREZFFiM1+AkRERI3EwUdERJbCwUdERJbCwUdERJbCwUdERJbCwUdERJbCwUdERJbCwUdERJbCwUdERJbCwUdERJbCwUdERJbCwUdERJbCwUdERJbCwUdERJbCwUdERJbCwUdERJZia/YTIH1RVQ2yokLTAEEAHJIIURSa/bSIysJ1TIfh4CPEUzKWYmkEN9KIJGUomrbzMUkQ0N3mQMDnxqDfjfYWRxOfKVFhXMdUKkHTdq0OspTVjTQeLm0glMhAEgW0OG1w2yVIu14ZK6qGdE5BKpuHomro9bowMehDn8/dxGdO9ArXMZWLg8+CMjkFU4sxPAsm4XHb0Oa0QRCK7wbSNA3JbB6JdA4nAh5cGPLDZZca8IyJ9uM6pkpx8FlMPCXjk+kwsjkVnR4HxBI2FK9TNQ1rmzKcdhHXT/VwtxE1HNcxVYODz0LiKRk3H4fgkAR4XPaqv95mJgdZ0fD+eC83GtQwXMdULZ7OYBGZnIJPpsM121gAgMdlh0MS8Ml0GJmcUpOvSXQYrmOqBQ4+i5hajCGbU2u2sdjmcdmRzamYWozV9OsSHYTrmGqBg88CVjfSeBZMotNTn904nR4HngY3EdzIlPU4VVUxPz+Pubm5ujwvMhe9rmMyHp7HZwEPlzbgcdsqOgCgFKIgwOu248FSHAFfoOjnq6qKhYUF/Ou//itWV1cxOTmJkZGRujw3Mg+9rWMyLg4+k4unZIQSGQR8rqKfe/d3nyGb3Xq1KwiAze6Ax+vD0bGTcLlbDn1sm9OG4EYG8ZRc8ACB3QMvFovB7/ejq6ur/N8UWU6169jn78DRsZNwOJyHPraUdUzGx8FnckuxNCRRKOn8pm3+zi443W7E19awHo1AURSMn3/j0McIggBJFLAUS+/bYBw08IaGhgAAsRibChVX8Tp2ubEWCSEaCgIAjo9PHvqYw9YxmQcHn8kFN9JocZb319zTN4iO7h7E/BFMP5xCOvWypMe1OG0IbqRxdsAH4PCBR1SOatZxq8eL2elHeLmZKOlxr69jMh8OPhNTVQ2RpIyO1vJeuYZXl7ARX8PG+hoAoKO7p6THue0SIkkZ+byCxcUXHHhUE5Wu463HqniZ3Bp4LW2ekh6zvY5VVeOFrU2Kg8/EZEWFoml7rllYithadOe/RVFAa5u3pMdJooDNZBJ//T//F+JrEfj9frjdbqRSqYKP2dzcxMbGBjye0jZKZD2yomE1KqH7+JGyHvfk23s7/+1t92P4+OmSHieJAhRt6+4OLpGXMjMjDj4Tq/SaPCfPnkdHdw82N+J4fO825p48grfdX/QAFwBwudy48r3v4dH9KaTTabS3t8PpLHxAgcvlQm9vL44dO1bZkyXTy+QUPM1tlv04f2cXVE3Dxvoakok40qmXsDtKf9fIa1qZFwefiVV71Heb1wdRlJDP55FJp0oafKIk4srld3DtnbcwMzODr776CslkEl1dXXC59h+Rp2ka/H4/ent7q3uyZFqZnALncq7sx203vpnvvkUkuILnM9OYfOPtkh9fp7MmSAc4+EzMIYmQBAGKWt7uzu3G93JzE/l8HqIoltRHFFWDJAj/dtNPCadPn8bY2NjOAIxGowUHIFEhla7jbUPDo4iGVvBycxPx9SjaOw4/hWb3OiZz4uAzMVHcuvlmUlbQVsYRcduNz2azweNrx9DwaNHznwAgnVPQ3ebYc0CA3W7nAKSqVLqOtzldbnT19iMSXMHywnzRwXfQOiZz4d0ZTO7b5Q08WIqjx1v/IRNOZDA52H7oYeC5XG7PLlBJkjAyMoL333+/7s+PjEtv65iMje/4TG7Q78bUixg0TSvr5N9yaZqGUCSChcQ80iutyOVyyGazyOVykGV555/t/89ms4hEIpidnUU+n+fgo4JSqRSSoQU8eRZB98UzdV/Hiqph0M87s5sZB5/Jtbc40Ot1IZnN1/yK9rsls3l47Rp++bP/C4fDgf7+fnR2dsLhcEAURUiSBEmS4HA44PF4IIoihoaG0NXVhQ8//LBuz4uMKxaLYXZ2FtFoFEeOHMHlM8cbso57vS5etcXkOPgsYGLQh18/CqHVWZ8L/KqahkQ6h//w1gT+3cB/wS9+8QsAwOLiInp6etDb2wu7ff/GKplMoru7mxeoph2apmF1dRVzc3PIZrMYGRnB+fPnYbPZ4N9IN2Qdvz3C68eaHQefBfT53DgeaMPzSArd3uIHqZRrbVPGiYAHAZ8Lvd4ziMVi+Pbbb3H06FGEQiE8ePAAfr8fgUAALS2vTomIxWL4/ve/X9ddV2QMuVwOi4uLmJubg9vtxtjYGHp7e/esjUauYzI3Dj6LuDDkx2o8g81Mrqa7ijYzOTjtIi4M+QFsXeT37bffxvr6OiKRCI4dO4aBgQFEIhE8efIEbrcbgUAALpcLLpeLJ65bXCqVwvz8/M7egTfeeAPt7e0FP79R65jMjUd1Wkg8JePm4xAcklCTjcZmJgdZ0fD+eO++JpJKpfDTn/4UgiDA53t10er19XUEg0GEw2HcuHEDN27cgCTxslBW83q/Gx4ehttd2gEl2+tYzcnwtTjgKONqLAc5bB2TOXHwWUw8JeOT6TCyORWdHkdFrUTVNKxtynDaRVw/1VNwY7G2toaPPvoIHR0de87ZUxQFs7OzuHz5MtLpNI4dO4Zjx44demkzMr6D+t3Q0BBstvJ3PIXjL/G/f/4FevuHcKyvo6p1HFpdwn985xQGezrK/hpkTBx8FpTJKZhajOFpcBNetx1tTltJnU3TNCSzeSTSOZwIeHBhyA+X/fB3a/Pz8/jHf/xHDAwM7GzggsEgTp06hatXryKZTGJ+fh7Ly8sIBAIYGRmB11vaRbHJGF7vd6Ojo/v6XTk0TcPt27dhd7mRaeureh3/7uOfYC0cxJ/8yZ/g7Nmzezo0mRMHn4UFNzJ4sBRHKJGBJApocdrgtkt7LgulqBrSOQWpbB6KqqHX68LkYHtZBwDcvXsXt27d2rk10dLSEn74wx+io+PVK2xZlrGwsIDnz5/D4/FgZGQE3d3dPPDFwF7vdyMjI4f2u1LNzs5idXUV77zzDkRRLLiORQGIRqPo7u4+dB1/+eWX+OlPf4qenh50dHTg/PnzHIAmx8FHiKdkLMXSCG6kEUnKUHYtCUnYulxUwOfGoN9dUQNRVRWffPIJZmdn4Xa7Dz13T1VVrKysYHZ2FqqqYmRkBIODg+yABlJNvytmfX0dt2/fxrVr1/Z9zdfXcV5R8PjxY5w5e/bQdTw3N4ePP/4YGxsbOH78ODY3NyEIAgegiXHw0R6qunUfMk3bujr91gWnq3/XJcsyfv7zn2NhYQE//OEP0d/fX/Qx0WgUc3NziMVi7IA6V8t+V4gsy/jss89w9uxZBAKBQz9XVTVk5Bz+369/jQ8++ODQdRyNRvHRRx/B4XAgGAxifHwcmqYhHA5zAJoUT2egPURRqMvNNx0OB27cuIGHDx+ir6+vpMd0dXWhq6trpwP+9re/ZQfUmVLOv6sFTdNw79499PX1FR16wNY6dtolOCShaIfeXkvd3d1IJBJ48eIFhoeH0d/fj3w+j6mpKdy7d48D0ET4jo8Mgx1QP+rV7wp5veuVQlEU3Lx5s6RL4v34xz+G0+mEzWbDo0eP0N/fj66uV1dwyefze94BXr58mevOwPiOjwzD4XDg+PHjGB0dxcrKCr777js8evSIHbCBXu93v/d7v1ezflfI+vo6ZmZmcO3atZKHXrl6e3uxsrICv9+PsbExTE9Po7W1def3ZrPZ0NnZidXVVaTT6bo8B2ocDj4yHFEUMTg4iMHBwZ0OOD09zQ5YJ4ddP7PeZFnG3bt3ce7cuboO2EAggLm5OQBAS0sLhoaGMDMzg/HxcUiShFgshnQ6jQ8++ABjY2N8t2dwHHxkaOyA9dOofldIuV2vGn6/H7urz3bve/78OdxuN3w+H/7wD/9wzyk4ZFwcfGQKbW1tmJiYwMmTJ7GwsICvvvqKHbBC5V4/s17m5uYgyzLeeOONun8vr9eL1w936Ovrw5dffomrV6/ij//4j6u+NBrpBwcfmQo7YOWa0e8KaUTX2621tRWSJEFRlD27Nv/qr/4KoVAIsixz8JkIBx+ZEjtgaZrZ7wppVNfbTRRFdHd3I5VKYXNzc8+uzRcvXuDOnTu4evUqXziZBAcfmR474H7N7neFNLLrvS4QCOCzzz7Dm2++iXfeeWfnHd6RI0ewtraGb7/9FufOnWvoc6L64Hl8ZDlWPh+w0efflauS8/UKKec8PgCIRCJIJBIYGRnZtxby+Tw+//xzHD9+HIODg1U9L2o+Dj6yLCtdF7Se18+slcOuw1mJcgdfMYlEArdu3cKVK1fQ1tZWk69JzcFdnWRZZu+Aeux3hTSj65XL6/Xi9OnT7H0moL+fAKImMFMH1Gu/K6SZXa9cR44cQTQaZe8zOA4+ol2MfD6gXs6/K1cjz9erhcnJSXz++edYWlpi7zMoDj6iAxjpfEA9nX9Xrkafr1cLNpsNly5dwq1bt9De3s7eZ0AcfESH0GsHNFK/K8QIXa8Q9j5jM85PCVGT6aEDGq3fFWKkrlcIe59xcfARlakZHdCo/a4Qo3W9Qtj7jImDj6hCjeiARu53hRix6xXC3mdMHHxEVap1BzRDvyvEyF2vEPY+4zH+TxKRjlTTAc3S7woxQ9crhL3PWDj4iOqgnA5otn5XiFm6XiHsfcbBwUdUR4d1wNbWVjx//txU/a4QM3W9Qtj7jIODj6gBtjvgwMAAHj16hJs3b2J9fR2XL1/G1atX0dra2uynWDdm7HqFsPcZAwcfUQO83u/+6I/+aOcd3+eff27Y64IWY+auVwh7n/5x8BHVUbF+Z9TrgpbK7F2vEPY+fePgI6qDWCyGubk5RCKRov3OSNcFLYcVul4h7H36xsFHVCMHnX937ty5ks+/0+t1QSthpa5XCHuffnHwEVUpn8/jxYsXNT3/Tg/XBa2UFbteIex9+sTBR1ShRpx/Z8T7A1q16xXC3qc/HHxEZSqn39WKUTqglbteIex9+sPBR1SCavtdrei5A7LrFcbepy8cfESHqEe/qxU9dUB2veLY+/SDg4/oAEa6fqYeOiC7XmnY+/SBg49ol2b0u1ppVgdk1ysde58+cPCR5eml39VKIzsgu1752Puaz5g/2UQ1oOd+Vyv17IDsepVj72suDj6yHCP1u1qpRwdk16sOe1/zcPCRZRi539VKrTogu1712Puah4OPTM1s/a5WqumA7Hq1w97XHNb+6SfTskK/q5VyOiC7Xu2x9zUeBx+ZihX7Xa2U0gHZ9eqDva+xOPjIFNjvaqdQB+zo6MDy8jLeffdddr0aY+9rLA4+MqzX+93w8DD7XQ3t7oArKyv4+7//e3R2duLFixdNvy6oGbH3NQ63EGQ47HeNpWkalpaW8O677+Lo0aNNvy6ombH3NQYHHxkG+11z7O56oig2/bqgZsfeV38cfKR77HfNU+h8PaPcH9CI2Pvqj4OPdInn3zVfKefr6fn+gEbG3ldf3IqQrrDf6UMl5+vp6f6AZsDeVz8cfKQL7Hf6Us35enq4P6BZsPfVBwcfNRX7nf7U6jqc7IDVY++rDw4+ajj2O/2qx3U42QGrw95Xe9zSUMOw3+lbI67DyQ5YGfa+2uLgo7pjvzOGRl6Hkx2wfOx9tcPBR3XDfmcczbq/Hjtg6dj7aoeDj2qK/c549HB/PXbA0rD31Qa3RlQT7HfGpMf767EDHo69r3ocfFQV9jtj0/P99dgBC2Pvqw4HH1WE/c74mtX1ysUOuB97X3U4+Khk7HfmoYeuVy52wL3Y+yrHLRYVxX5nLnrseuViB9zC3lcZDj4qiP3OnPTc9cp1WAfs6Oho9tNrCPa+8nHw0T7sd+ZllK5XroM6YC6XQzgchqIopt4NyN5XPkHTNK3ZT4Ka76B+NzQ0xH5nIrIs47PPPsPZs2cNu4uzHKFQCD/5yU8wPj5uiQ64nSPY+4rj4LO4XC6HxcXFnX43OjrKfmdCmqbh9u3baG1txZkzZ5r9dBpCURTcvHkT165dw/z8PJaXl03fAe/evQtJktj7iuDLeYtiv7MWM3W9clnpfED2vtJw8FlMLBbD7OwsotEo+51FxGIxU3a9clnhfED2vtJw8FnAQf3u/Pnz7HcWkMvl8M033xjqfL16M/v5gDy/rzhu+Uzs9X7H8++sZ2pqytDn69WbWc8H5Pl9h+PgMyH2OwKA2dlZy3a9cpmxA7L3FcbBZyLsd7Rtey1cvXrV0l2vXGbqgOx9hXHwGRz7Hb1uu+tNTk7yhU+FzNIB2fsOxq2jQbHfUSHserVl9A7I3rcfB5/BsN/RYdj16sfIHZC9by8OPoNgv6Ni2PUaw4gdkL1vLw4+HWO/o1Kx6zVeoQ549OhRDA8P664Dsve9wi2oDrHfUbnY9ZrLKB3wyJEjWFtbs3zv4+DTEfY7qgS7nn4c1AHb2towOjqqmw44MTFh+d7HwacD7HdUKXY9fdJzB2Tv4+BrGvY7qha7nv7ptQNavfdxK9tg7HdUK+x6xqK3Dmjl3sfB1yDsd1RL7HrGpacOaNXex8FXZ+x3VGvseuaghw5o1d7HwVcH7HdUL+x65tPsDmjF3sctcQ2x31G9seuZW7M6oNV6Hwffv1FVFblcrqJXV+x31AjsetbRjA5opd7HwYetXZNffvkl4vE4fvCDH5T8OPY7ahR2PWtqZAe0Uu/j4APw6NEj3L9/H6IoIhqNoqurq+Dnst9Ro7HrUaM6oFV6n6G21qqqQVZUaBogCIBDEiGK1b3lX1pawqeffor+/n7E43E8ePAA169f3/d57HfULOx6tFu9O2A5va8e2+RG0P3gi6dkLMXSCG6kEUnKUDRt52OSIKC7zYGAz41BvxvtLY7yvnY8jn/6p39CV1cX7HY7Ojs78eTJE7z55ps7b/PZ76iZ2PWokHp2wMN6Xz23yY0iaNquZ60jqxtpPFzaQCiRgSQKaHHa4LZLkHa9mlBUDemcglQ2D0XV0Ot1YWLQhz5f8d1BmUwG//AP/4BcLge/37/z6ysrK7h06RJGR0f39Lvh4WHuZqKGisViuH37Nq5evcq1VwFFUXDz5k18+OGHzX4qDaGqKlZWVjA7OwtVVavugIlEArdu3cKVK1fQ1tZW921yI+lu8GVyCqYWY3gWTMLjtqHNaSvplYumaUhm80ikczgR8ODCkB8u+8F/4ds/EEtLS+jr69vzNUKhEGZmZvDuu+/i5MmTGBoaYr+jhsvlcvj0009x9uxZ7uKskNUG325ra2uYnZ1FLBarqgMuLi4iGosj09ZX121yo+lqix5PyfhkOoxsTkWg3QWxjLfqgiDA47Kj1WnD80gKq/EMrp/qOfCt9tdff43nz59jaGgIAJDP5xGNRhEMBuFwONDR0bGzWIiagV2PqtHZ2YnOzk68fPkSc3NzFXdAT2cvbkcEZCOpum6TG003x0XHUzJuPg5BEIBur7OsP+DdREFAt9cJQQBuPg4hnpL3fPzJkye4c+cO+vv7kc1m8eLFC9y/fx8vX77E2NgYxsfHcezYMXzzzTdQVbUWvzWismx3vdOnTzf7qZDBtba2YmJiAtevX0drayu++uor3Lp1C+FwGMV29jVqm9wMuhh8mZyCT6bDcEhbrxBqweOywyEJ+GQ6jExOAQCEQiH88z//M7xeL+bm5vDtt98CAM6ePYvR0dGdA1paWlqQSCSwvLxck+dCVKrt8/UuXbrE8/WoZrbPB3zvvfcwNDSE7777Dv/yL/+ChYUFKIqy7/MbtU1uFl3s6pxajCGbU9Htre016TwuOyKJLKYWYzjT5cCPf/xjbGxswG6377ztLxR+vV4v7t69u7M7lKjeeL4e1dvu8wG3O+BB5wM2Ypv89kjh86XrremDb3UjjWfBJALtrrp8/U6PA0+Dm/jyV7/B8uIi+vr64PP5oCgKQqEQNE2Dpmk7sVYURUiSBFVV8fDhQ5w+fRonTpyoy3Mj2o1djxqpUAfsGBxtyDZ5uLMNAV99vkcxTR98D5c24HHbKt5/XIwoCPC67ej+3h/gv/ynH0FRFOTz+X3/3v4nmUxiZmYG9+7dQzweRzwer8vzItqN5+tRs2x3wJMnT2J5eRkPl+MN2SY/WIoj4GvOi7ymDr54SkYokSl56ifiMay8mMfmRhyqqsDucMHf1YWjoycP7SFtThuCGxnkBDvava0Hfo4sy3jy5Anu37+/swGKRqPo7Oys6PdGVCpeh5P0wOFwwN87gPXgSsnb5Lu/+wzZbAYAcP7NK3C3Hrx9fd32NjmekptylGdTB99SLA1JFEo6JyQaWsXMdw+haUBLWxvaPD5ksxmElhcxdGzs0A2GIAiQRAFLsfS+P+Ttgff1119DlmV0d3fX/f5XRNvY9UhPytkmJ+LrO0MPACKhFRwZOV7S9zlsm9wITR18wY00WpzFn4KiKJh/+h00Deju7cPo6bM7fzHp1EuIJVyZoMVpQ3AjjbMDPgAceKQP7HqkJ6VukwEgElwFALR6PHi5uYlocLXkwQfs3yY3UtMGn6pqiCRldLQWn/abGzHk83kAwMDRkT2vRtwtpb21dtslRJIyMpksnj17yoFHTceuR3pSzjZZVVWsR0IAgKOjJ/D02/vIZjNIxNfhbe8o6fttb5NVVWv4ha2bNvhkRYWiaXuu81ZIPpfb+W9nhbuDJFHA4tIy/sfP/g+UXBadnZ2w2+1YWFgo+JhwOIz79+9jbW2tou9JVEgikcDjx49x/vx53Llzp9lPx5QURcGTJ08Ovc0YvaIKEhStq6Rtcmwtgnw+D7vDAW97B/yd3YiEVhEJrpY8+CRRgKJt3d3BJTb2UmZNG3zlXCHUZn91AmU2nS45oL7O7/fj+Ftv49n0I4iiiM7OzkMv4KooCgYGBnD06NGKvh/RQXK5HBYWFnD9+nV0d3c3++mYlqIomJub489viXKagJmFdEmfGw2uAAD8nVt3gejo7kUktIr1SAjDJ06XdZBWM64W3bTBV86Rsh6fHzabDfl8HssLc3saXzaTht3hLOkP2t3ixh987wPceO/38fDhQ9y/fx+CIKCnp+fAAfjy5Ut0dXWht7e39CdLVMTXX3+NM2fOYHx8vNlPxdQURYHf7+fPb4kyOQVYWCr6eflcDvH1KAAgvLqM8OqrK1zl83nE1iLo7C79z7wZtzRt2uBzSCIkQYCiFt/dKUkShk+cxsx3DxEJreLly020eXzIyVlsxNZw6Z3fLzr4FFWDJAhbN0psbcVbb72FiYmJkgYgUa1sd71Tp041+6kQ7VHqNnktHISqapAkCT7/q92aqZcvkUmnEA2ulDT4dm+TG61pg08Ut25YmJQVtJVwFFFXbx8cTieWF54jmYgjmnoJu8OFnr7Bko7qTOcUdLc59kTUVg5AaiCer0d6Vuo2ORLaOpqzt38QR8dO7vx6Ir6OR1N3EF+PIifLsDsOP0jmoG1yozT1dIaAz40HS/GSBh8AeNs7Sg6nr0tl8xjrbjvwY4UG4EEXbyWqxPb5eufOneP5eqRbpWyTz15888Bf97Z34O13b5T8vQ7bJtdbU192DvrdUFSt6O0xqqVpGhRVw6D/8A3O9gD8y7/8S0xOTkIUxZJO5CQqZmpqCv39/exNpEv5fB5zc3N4/uguUpmsbrbJ9dLUd3ztLQ70el1IZvM1u/XFQZLZPHq9rpKvELA9ACcnJ2G31+95kTWw65FepdNpzM/P48WLF+ju7saVy+dxeyWru21yrTU9NEwM+rCZzkOt0ysMVdOQSOcwOdhe9mNbWlo4+KgqvL8e6VE8Hsc333yDTz/9FJqm4dq1a7h06RL8fr+ut8m10vS7M/T53DgeaMPzSKrm934CgLVNGScCnqbd/oKsi12P9ETTNASDQczNzSGTyWB4eBjnzp2DzbZ3DFhhm9z0wQcAF4b8WI1nsJnJ1fTt9WYmB6ddxIUhf82+JlGp2PVID/L5PF68eIH5+Xk4nU6Mjo4iEAgcevyC2bfJuhh8LruE66d6cPNxqGZ/0JuZHGRFw/vjvXDZeWoCNRa7HjXb6/3u4sWL8PtLGzjb2+RfPVxBKp1Gr99b9fPR0zZZF4MP2DrQ5f3xXnwyHUYkkUWnx1HRjRBVTcPapgynXcT7471Ni6dkXTxfj5opHo9jdnYWkUgEQ0NDuHbtGlpaWsr+Ou0tDvTmQ3gcFyBKTlNtk3Uz+ICtP+gPJ/owtRjD0+AmvG472py2kk4p0DQNyWweiXQOJwIeXBjyN/1VBVkPux41Q6n9rhxzc3OwqzL+879/C/eWN0y1TdbV4AO23mK/PdKF4c42PFiKI7iRgSQKaHHa4LZLey6lo6ga0jkFqWweiqqh1+vC2yNdPJCFmoZdjxqpkn5XilgshmfPnuHq1atwO+2m2ybrbvBtC/hcCPgCiKdkLMXSCG6kEUnKUHYdYisJW5fYGetuw6DfrYu30GRd7HrUKNX0u2J277XYvYvUTNtk3Q6+be0tDrS3OHB2wAdV3bp3k6ZtXdHbIYlNuc4b0evY9agRatXvDnPv3j309fUhEAgc+HEzbJN1P/h2E0Wh4TcsJCqGXY/qqR79rpDt73Hp0qWSPt+o22RDDT4iPWLXo3qoV78rZHfXM/teCw4+oiqw61Gt1bPfFVKo65kVBx9Rhdj1qJYa0e8KKdb1zIaDj6gCsiyz61HVGtnvCim365kBBx9RmTRNw71799j1qGKN7neFWKnr7cbBR1Smubk5dj2qSDP6XSFW63q7cfARlWF9fZ1dj8rWzH5XiNW63m4cfEQlkmUZd+/eZdejkuih3xVixa63W/P/BogMgF2PSqWXfleIVbvebhx8RCVg16Ni9NTvCrFy19uNg4+oCHY9Oowe+10hVu56u3HwER2CXY8Ooud+V4jVu95u+v1bImoydj16nd77XSHsentx8BEVwK5H24zQ7wph19uPg4/oAOx6BBir3xXCrrcfBx/Ra9j1rM2I/a4Qdr2DGe9vkqiO2PWsy6j9rhB2vcI4+Ih2YdezHiP3u0LY9Q7HwUf0b9j1rMUM/a4Qdr3DcfARgV3PKszU7wph1yvOPH/bRBVi1zM/s/W7Qtj1SsPBR5bHrmde2WwWjx8/NlW/K4Rdr3QcfGRp7HrmFI/H8fTpUzx8+BAnTpwwVb8rhF2vdBx8ZFnseubyer87cuQILly4gDNnzjT7qdUdu155OPjIktj1zKNQv1NVFc+ePWv206s7dr3ycfCRJbHrGZ8Zz78rF7teZTj4yHLY9YzNzOfflYtdrzIcfGQp7HrGZIXz78rFrlc5664ashx2PeOxyvl35WLXqw4HH1kGu55xsN8Vxq5XPQ4+sgR2PWNgvyuOXa96HHxkeux6+sZ+Vzp2vdrgyiJTY9fTL/a78rDr1Q4HH5kau57+sN+Vj12vtjj4yLTY9fSF/a5y7Hq1xcFHpsSupw/sd9Vj16s9rj4yHXa95mO/qw12vfrg4CPTYddrHva72mHXqx8OPjIVdr3mYL+rPXa9+uHgI9Ng12ss9rv6YderL65QMgV2vcZhv6svdr364+AjU2DXqz/2u/pj12sMDj4yPHa9+mK/axx2vcbg4CNDY9erD/a7xmPXaxyuYjIsdr3aY79rDna9xuLgI8Ni16sd9rvmYddrPA4+MiR2vdpgv2s+dr3G4+Ajw2HXqw77nX6w6zUHVzoZCrte5djv9IVdr3k4+MhQ2PXKx36nP+x6zcXBR4bBrlce9jv9YtdrLg4+MgR2vdKw3+kfu17z8aeBdI9drzj2O2Ng19MHDj7SPXa9wtjvjINdTz84+EjX2PUOxn5nPOx6+sHBR7rFrrcX+51xsevpC39iSJfY9V5hvzM2dj394eAjXWLXY78zA3Y9feLgI92xetdjvzMPdj194uAjXbFq12O/Mx92Pf3iTxXphhW7HvudObHr6RsHH+mGlboe+515sevpHwcf6YJVuh77nfmx6+kfBx81ndm7HvuddbDrGQN/8qipzNz12O+shV3PODj4qKnM2PXY76yHXc9YOPioaczW9djvrItdz1g4+KgpzNL12O+IXc94+NNJDWeGrre737lcLoyMjLDfWRC7njFx8FHDGbnrsd/RNnY94+Lgo4Yyatdjv6PXsesZFwcfNYzRuh77HRXCrmds/AmmhjBS12O/o8Ow6xkfBx81hBG6HvsdFcOuZw4cfFR3eu967HfG9/HHHyMaje75NUVR8OjRIwSDwX2fPzAwgPfff7/s78OuZw4cfFQ1VVULDjS9dj32O3Pp6enBysrKnoGkKAo8Hg86Ojr2fO7Kygr6+vr2fQ1N0wCg4C5tdj3z4E85VUVRFHz00Uc4d+4cTp48uWejoceux35nThMTE7h//z40Tdvz4kUURUiStPP/sizD6XTi5MmT+77G06dPMT09jffeew9tbW17PsauZy78G6SqbG5uYm1tDb/5zW/w29/+FrIs73xMT10vnU7j8ePH+M1vfoNYLIaLFy/iypUr6Ovr49AzgdbWVpw7dw6RSOTQzwuHw7h8+TKcTue+j62srGB+fh5/+7d/i8XFxZ1fZ9czH77jo6okEgmIoojBwUE8e/YMoVAIN27cgCAIuuh67HfWsf2uL5/PH7jLWpZlOByOA9/tAcDq6ioGBgagqio+/vhjvPnmm7h48SK7nglx8FFV1tbWIEkSBEFAf38/4vE4/uZv/gYtLS348MMPm9L12O+saftd34MHDw5seOFwGFeuXDnw3V4ul0MikUB/fz8EQcDAwAC+/vpr3L17F8eOHWPXMxluCagqq6ure95B+Xw+LC8vIxgM4vHjx/D7/QduaOqB/Y52v+vb/fde7N1eIpEA8OrAFkmS0NHRgdu3b0NVVSwvL2NoaKj+vwFqCA4+qko4HIbX6935/2AwCEEQ8MYbb+zZ9dnZ2Vm358Dz72jb7nd9PT09O79+2Ls9YGvwbR/VCWy9iJqZmcH58+fhdDr37PrcfbAMGRMPbqGKpdNpZDIZ2O12AFsHugSDQYyOjkKSJPT39yOfz+Pv/u7vMD09vWfDUgvxeBzffPMNPv30U2iahmvXruHSpUscehY3MTEBTdOQz+cBFH+3BwCRSGRnHQPA/Pw8Ojo64Pf70dLSgoGBAdy+fRu/+MUvkEwm6/57oPriOz6q2PbuIWCrkczOzmJ4eHjPq+r29nYIgoBf/vKXsNlsGBsbq+p7st9RMdvv+qampgAUf7cHbO2p2N5lHwwGIcsyRkdHdz4uSRIGBgbw7NkzyLKMP/3TP+VpDQbGrQVVbGNjA8DWMNp+hdze3r7z8VQqhbW1NXg8HvzgBz/A8PAwMpkMXC5X2d+L/Y7KMTExgampKWSz2aLv9jRNQzgcRldXF5LJJFZWVnDmzJmdwaaqKqLRKGRZxsWLF3HhwgUOPYPj4KOKhUIhuFwuBINB5PN5DA4OAtg78L7//e9jeHgYNpsNMzMz+N3vfoc/+7M/K3n4sd9RJbbf9X311Vf40Y9+dOi7vZcvXyKfz0PTNMzMzOzstdg98E6ePIkLFy5w7ZkEBx9VbHV1FYqiIBgMYnx8HJlM5sCBB2w1lF//+tc7u0TPnDlz6NeOxWKYm5vj+XdUsYmJCYyNjR36bg94tct+e6+Fz+dDOBzmwDMxDj7aQ1U1yIoKTQMEAXBIIkRx/+5ERVEQDoextraGvr4+hMPhAwcesPWK+pe//CV8Ph/sdjvu3LmDU6dO7Ts6jv2OakVVNdhdLfiv/+2/QxMFqKp24DoGtg6SikajEEURHo8HKysrHHgmxy0KIZ6SsRRLI7iRRiQpQ9l19KUkCOhucyDgc2PQ70Z7iwPA1qvk58+fw+12w+1246233to38ICtg15u3rwJRVF2LhYciUSwuLiIY8eOAWC/o9qoZB0DwLNnz/DixQuMj49jdHSUA88COPgsbHUjjYdLGwglMpBEAS1OGzpaHZB2vTJWVA1JWcGDpTimXsTQ63VhYtCHteVltLW14c///M8xOjp64LsyTdPwxRdfIBwOY2BgYOfXfT4f7ty5g56eHjx//pz9jqpSzTru87kxMzODGzdu4L333uP6swhBq/XJVaR7mZyCqcUYngWT8LhtaHPaSnp3pWkaktk8EukcTgQ8OD/gg9tpL/j59+7dwxdffIEjR47s+frJZBL37t3D+fPnMTk5ieHhYfY7Klut1vG5fi9aXI6ijyPz4OCzmHhKxifTYWRzKjo9DogV7E5UNQ1rmzKcdhHXT/Xs2W20bWFhAT//+c8xMDAAm80GTdMQi8UQDAaRy+XgcDgwOTmJDz74oBa/LbKYRq1jMicOPguJp2TcfByCQxLgcRV+p1aqzUwOsqLh/fHePRuN9fV1fPTRR2hvb4fD4UAkEkEoFILdbkcgEIDf74emaVhZWcFf/MVf7LnkGVExjVrHZF48C9MiMjkFn0yHa7axAACPyw6HJOCT6TAyOQXA1nl3v/rVryBJEiKRCO7du4dkMonR0VGMj4+jo6MDgiDs3CD08ePHNXkuZA2NWsdkbhx8FjG1GEM2p9ZsY7HN47Ijm1MxtRiDoij42c9+hunp6Z0beZ49exZjY2P77mgNAN3d3Xj48CEymUxNnxOZVyPWMZkfj+q0gNWNNJ4Fkwi0l3+psFJ0ehx4GtzE8nf3cOvWLfT09KCjowOSJGFtbW3Pxak1TYMkSTvv+BKJREkntBM1ah0Pd7Yh4KvP9yB94OCzgIdLG/C4bRUdAFAKURDgddvh9o7jr//6Laiqinw+D0VR9v07l8tBluWdf7LZLG/zQiVp1Dp+sBRHwMe7rZsZB5/JxVMyQolM0Vewd3/3GbLZrV2OoijCbnegzedD3+BReHztRb9Pm9OG4EYGedGB9jYeIEC1Veo6Bvau5d0mL7+N1jbPoY/dXsfxlMwDXUyMg8/klmJpSKJQ8lVQ/J1dsNkd2NyIYy0cwnokhLHxSXT1HP4KWBAESKKApViaGwyquXLXMbC1ll3uV+eH7r7fXiFcx9bAwWdywY00Wpyl/zX39A2io7tn60r13z1ENBTE/JPH8Hd2F90l2eK0IbiRxtkBX7VPm2iPctcx8Gotl4vr2Pw4+ExMVTVEkjI6Wst/5SoIAgaPjiIa2rrl0OZGDO0dXYc+xm2XEEnKh14QmKhcla7j8OoSEvH1nf8/dvxUSY/jOjY/Dj4TkxUViqbtuWZhORy77pmXz+WKfr4kClC0rbs7uEQesEK1Uek6jq1F9/x/qYOP69j8OPhMrNpr8si7zq+zldBHavV9iXardD2dPHu+ol2d1X5f0j+ewG5i1Rz1rWkalhZmAQA2mw0eX+lXrefdhKiWmrWeuI7Ni+/4TMwhiZAEAYpa+m6i8OoS1qNhbG7EkUmnIAjA8Mnxks61U1QNkiDAIfH1FNVOJesY2N/4uvsGip7OAHAdWwEHn4mJ4tbNN5OygrYSj4iLrUV3zuPr7OlF/9AxtHlLO7otnVPQ3ebgAQFUU5WsY2B/4/O2d5Q0+LiOzY+Dz+QCPjceLMWLbjAuvnOt6u+VyuYx1r3/mpxE1Sp1HQPVr2WuY/Pje3mTG/S7oaga6n33KU3ToKgaBv3uun4fsiauY6olDj6Ta29xoNfrQjKbr+v3SWbz6PW6eLULqguuY6olDj4LmBj0YTOdh1qnV8uqpiGRzmFysL0uX58I4Dqm2uHgs4A+nxvHA21Y25Tr8vXXNmWcCHh4KxeqK65jqhUOPou4MOSH0y5iM1P8Cizl2Mzk4LSLuDBU+nl+RJXiOqZa4OCzCJddwvVTPZAVrWYbjc1MDrKi4fqpHrjsvLQT1R/XMdUCB5+FtLc48P54LzQNiCSyFbcSVdMQSWShacD74708EIAaiuuYqiVo9T4+mHQnk1MwtRjD0+AmvG472py2ku5zpmkaktk8EukcTgQ8uDDk5ytkahquY6oUB5+FBTcyeLAURyiRgSQKaHHa4LZLey4Lpaga0jkFqWweiqqh1+vC5GA7DwAg3eA6pnJx8BHiKRlLsTSCG2lEkjKUXUtCErYuFxXwuTHod3N3EOkW1zGVioOP9lDVrfuQadrW1ekdkshrFpLhcB3TYTj4iIjIUnhUJxERWQoHHxERWQoHHxERWQoHHxERWQoHHxERWQoHHxERWQoHHxERWQoHHxERWQoHHxERWQoHHxERWQoHHxERWQoHHxERWQoHHxERWQoHHxERWQoHHxERWcr/B+Kxj21+0GQVAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# feature R correlates with other features\n",
    "model.add_edges_from([(\"R\", \"B\"), (\"R\", \"C\"), (\"R\", \"D\"), (\"R\", \"E\")])\n",
    "nx.draw_circular(\n",
    "    model, with_labels=True, arrowsize=30, node_size=800, alpha=0.3, font_weight=\"bold\"\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Then, parameterize our graph to create a Bayesian network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pgmpy.factors.discrete import TabularCPD\n",
    "\n",
    "# add CPD to each edge\n",
    "cpd_a = TabularCPD(\"A\", 2, [[0.7], [0.3]])\n",
    "cpd_r = TabularCPD(\n",
    "    \"R\", 3, [[0.6, 0.2], [0.3, 0.5], [0.1, 0.3]], evidence=[\"A\"], evidence_card=[2]\n",
    ")\n",
    "cpd_b = TabularCPD(\n",
    "    \"B\",\n",
    "    3,\n",
    "    [\n",
    "        [0.1, 0.1, 0.2, 0.2, 0.7, 0.1],\n",
    "        [0.1, 0.3, 0.1, 0.2, 0.1, 0.2],\n",
    "        [0.8, 0.6, 0.7, 0.6, 0.2, 0.7],\n",
    "    ],\n",
    "    evidence=[\"A\", \"R\"],\n",
    "    evidence_card=[2, 3],\n",
    ")\n",
    "cpd_c = TabularCPD(\n",
    "    \"C\",\n",
    "    2,\n",
    "    [[0.7, 0.2, 0.2, 0.5, 0.1, 0.3], [0.3, 0.8, 0.8, 0.5, 0.9, 0.7]],\n",
    "    evidence=[\"A\", \"R\"],\n",
    "    evidence_card=[2, 3],\n",
    ")\n",
    "cpd_d = TabularCPD(\n",
    "    \"D\",\n",
    "    3,\n",
    "    [\n",
    "        [0.3, 0.8, 0.2, 0.8, 0.4, 0.7],\n",
    "        [0.4, 0.1, 0.4, 0.1, 0.1, 0.1],\n",
    "        [0.3, 0.1, 0.4, 0.1, 0.5, 0.2],\n",
    "    ],\n",
    "    evidence=[\"A\", \"R\"],\n",
    "    evidence_card=[2, 3],\n",
    ")\n",
    "cpd_e = TabularCPD(\n",
    "    \"E\",\n",
    "    2,\n",
    "    [[0.5, 0.6, 0.6, 0.5, 0.5, 0.4], [0.5, 0.4, 0.4, 0.5, 0.5, 0.6]],\n",
    "    evidence=[\"A\", \"R\"],\n",
    "    evidence_card=[2, 3],\n",
    ")\n",
    "model.add_cpds(cpd_a, cpd_r, cpd_b, cpd_c, cpd_d, cpd_e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next, generate sample data from our Bayesian network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Generating for node: B: 100%|██████████| 6/6 [00:00<00:00, 192.64it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      A  R  B  C  D  E\n",
      "0     0  1  2  1  0  0\n",
      "1     1  2  2  1  2  0\n",
      "2     0  0  2  0  2  1\n",
      "3     1  2  2  1  0  1\n",
      "4     1  2  2  1  0  1\n",
      "...  .. .. .. .. .. ..\n",
      "9995  0  0  2  0  1  1\n",
      "9996  1  2  1  0  0  1\n",
      "9997  0  0  2  0  1  0\n",
      "9998  1  0  2  1  0  1\n",
      "9999  0  0  2  0  2  1\n",
      "\n",
      "[10000 rows x 6 columns]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "from pgmpy.sampling import BayesianModelSampling\n",
    "\n",
    "# sample data from BN\n",
    "inference = BayesianModelSampling(model)\n",
    "df_data = inference.forward_sample(size=10000)\n",
    "print(df_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now we are ready to learn the TAN structure from sample data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Building tree: 100%|██████████| 15/15.0 [00:00<00:00, 5215.93it/s]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAwzklEQVR4nO3d6XMbd37n8U934yYBECRBghQoHhJFnZRk2R55KGlmNDO2RhtXUpOqGacSVzapPMn/kN1K7YPUVu1uKk93n2SrcjzI2F7HqplE5Snb0cyItseWZF2mDlIiBZIgQRIgAOLu7n3AIU1KBHEQRx+fV9XU6CDAltXqL7rf+DUEVVVVEBERmYTY7A0gIiJqJA4+IiIyFQ4+IiIyFQ4+IiIyFQ4+IiIyFQ4+IiIyFQ4+IiIyFQ4+IiIyFQ4+IiIyFQ4+IiIyFQ4+IiIyFQ4+IiIyFQ4+IiIyFQ4+IiIyFQ4+IiIyFQ4+IiIyFUuzN4C0RVFU5GQFqgoIAmCTRIii0OzNIqoI92PaDQcfIZbKIRRNI7yaRiSZg6yqm78nCQL8rTYEvE4EfU60uWxN3FKi4rgfU7kEVd2yd5CpzK+mcSe0ioV4BpIowGW3wGmVIG15ZSwrKtJ5GalsAbKiotvjwImgFz1eZxO3nOgb3I+pUhx8JpTJy7j5LIpH4STcTgta7RYIQunLQKqqIpktIJ7O41DAjdN9PjisUgO2mOhF3I+pWhx8JhNL5fDRxCKyeQUdbhvEMg4Uz1NUFcuJHOxWERcPd/GyETUc92PaCw4+E4mlcrh6fwE2SYDbYd3z8yUyeeRkFW8c7eZBgxqG+zHtFZczmEQmL+OjicWaHSwAwO2wwiYJ+GhiEZm8XJPnJNoN92OqBQ4+k7j5LIpsXqnZwWKD22FFNq/g5rNoTZ+XaCfcj6kWOPhMYH41jUfhJDrc9bmM0+G24WE4gfBqpi7PTwRwP6ba4eAzgTuhVbidlqreAFAOURDgcVpxOxSry/MTAdyPqXa4gN3gYqkcFuIZBLyOkl974/o1ZLPrr3YFAbBYbXB7vOg/OAKH07XrY1vtFoRXM4ilcnyDANVcJfvxhngsirmZJ0isxqAoMqw2B3ydneg/MAJR3Pk1P/djc+AZn8GFomlIolDW+qYNvo5OdO/rgyRZsLIUwdSD+yUfIwgCJFFAKJrey+YS7ajS/XhpYR73b/0W0eUl2BwOdHb3wuFyYWH2GRS5+BtYuB+bA8/4DC68mobLXtlfc1dPEO3+LkR9EUzcuYl0aq2sx7nsFoRX0zi+z1vNphIVVcl+LMsynjz8GqoK+Lt7cODI8c2BmU6tQZR2X6zO/dj4OPgMTFFURJI5tLdUdslmcT6E1dgyVleWAQDt/q6yHue0Sogkc1AUlTcEppqpdD9OrEZRKBQAAPv6h7adJTpdLSUfz/3Y+Dj4DCwnK5BVdds9C8sRXV7a/LEoCmhp9ZT1OEkUIKvrd8V3iLwFFNVGpftxIZ/f/LHdWfm9OLkfGx8Hn4FVe0+ekeOn0O7vQmI1hvu3foupB/fgafOVfIPLXr8v0U4q3Z8s1m/W+GXTaThbSp/l1eL7kn7wzS0Gttd3fbd6vBBFCaoKZNKphn1foq0q3Z/cXh8slvXX9LPTU9h6V8ZsJg1FUeryfUk/eMZnYDZJhCQIkJXKLnduNL61RAKFQgGiKMLV6i75OFlRIQkCbBJfT1HtVLofS5KEwUNH8PjrO4gszGNtLYFWtxf5XBar0WWc+fZ3iy5nALgfmwEHn4GJ4vqHbyZzMloreGfnRuOzWCxwe9vQN3gANpu95OPSeRn+VhvfEEA1Vc1+3NndA5vdjtnpp0jGY1hKrcFqc6CrJ1jyXZ3cj42Pg8/gAl4nbodiZR0wXvr2hT19r1S2gIP+1rK/XpZlFAoFSJK0eWmKaCeV7McbPG3t8LS1V/y9Kt2PSX94tDG4oM+JmzNRqKpa0SL2SqmqioVIBNPxJ0jPtSCfzyObzSKfzyOXy23+b+PnhUIBqqpClmWcOnUKY2Njdds20r9G7seyoiLo4yezGxkHn8G1uWzo9jiQzBZqfkf7rZLZAjxWFb94/2ew2Wzo7e1FR0cHbDYbRFGEJEmQJAk2mw1utxuiKEIQBDx79gyDg4N12y4yhkbux90eB29XZnCstyZwIuhFIl2AUqf3Zyuqing6j/909gT+8i//Et3d3QCAZ8+eIZFIwOFwwOVyweFwwGq1QpIkCIKAZDKJ9vZ29PT01GW7yFgatR+PBtvq8vykHRx8JtDjdWI40IrlRK4uz7+cyOFQwI2A14Fjx45hbGwMra2tOHLkCAqFAm7fvo2pqSmkUtuXRESjUZw5c6aul67IOBq5H5OxcfCZxOk+H+xWEYlMvvQXVyCRycNuFXG6zwdg/Sa/r732Gnp7exGPxzEwMIDR0VE4HA48ePAAExMTiMViyGQycDgcGBgYqOn2kLE1aj8mY+PgMwmHVcLFw13IyWrNDhqJTB45WcXFw11wWL95i7jFYsEPf/hD2O12rK6uwmq1ore3FydPnkRnZydCoRCuX7+Onp6eXddTET2vkfsxGRePOibS5rLhjaPdUFUgEs9W3UoUVcXT8ApCoVm8fqRrxzcCuFwuXL58GWtra8hk1j/jTxRFdHZ24siRIwgGg2hpacEvf/lLPHjwANlsdk9/NjKPWu7HkXgW09MzGNvfwje0mIigqrwjndlk8jJuPoviYTgBj9OKVrulrM6mqiqS2QLi6Txa5QSu/vP/xoWx1/Dd734X+/bt2/E5njx5gp///OfYt2/f5lq9cDiMw4cP4/z580gmk3jy5AlmZ2cRCAQwNDQEj6e8m2KTudViPz4UcOP6B/+M5cUw/uAP/gDHjx+Hy1XePWlJvzj4TCy8msHtUAwL8QwkUYDLboHTKm27LZSsqEjnZaSyBciKim6PA6PBNrgtMv7u7/4O6XQa3d3d6O3txbe+9a0dB+CNGzcwPj6Ovr4+AEAoFMJbb72F9vZvFhfncjlMT0/j6dOncLvdGBoagt/v5xtfqKS97McBrwO/+c1v8O6776Krqwvt7e04deoUB6DBcfARYqkcQtE0wqtpRJI5yFt2CUlYv11UwOtE0OfcvBykqir+/u//HtFoFFarFR0dHYjFYggEAi8MQEVR8NFHH2FychJOpxOdnZ24fPnyjtuiKArm5uYwOTkJRVEwNDSEYDAIqcRtpohe2I8VBeFwGIGenqL7MQBMTU3hgw8+wOrqKoaHh5FIJCAIAgeggXHw0TaKsv45ZKq6fnd6myQWvWfhlStXEI1GMTMzg2AwiI6ODsTj8R0HYC6Xw5UrVzA9PY233noLvb29JbdlaWkJU1NTiEajGBgYwMDAAOz20vcMJVIUFVPTM/hff/u3+C9/9VcI+DuL7sdLS0t45513YLPZEA6HcfToUaiqisXFRQ5Ag+KbW2gbURTgsEpw2iQ4rNKuN+rt7e1FNpvF8PAwpqenkU6n4fF4sH//fmQyGfzrv/4r3nvvPYRCIVitVrz++us4e/Zs2QvWOzs78eqrr2JsbAy5XA4ff/wxbt26hXg8Xqs/LhmUIAC3b36JQiaF+3e+2nU/3mjKfr8fLpcLMzMzsFgs6O3thd/vx82bN/GP//iP+Pzzz19Yi0r6xDM+qtrTp0/xb//2bwgGg1hcXMTi4iKOHDmy7bLk1jPA733ve9u6XqXYAalcc3NzePfddxGJRNDd3Y233noLHR0dRb/+n/7pn2C322GxWHDv3j309vais7Nz8/cLhcK2M8BXXnmF+52O8YyPqubxeDb/8Xd1dcHpdGJmZuaFr9lYz1fuB4AWY7PZMDw8jO9///sIBoP4+uuv8cknn2B6ehqyLO/puck4VFXF559/DrfbDUEQ4HQ6cfPmzV0f093djVQqBUmScPDgQczMzCCdTm/+vsViQUdHB7LZ7LZfJ33i4KOqbVwi2rhoMDAwgHg8juXlZQDrr5JnZmbQ09ODn/70p9teQe+FKIoIBoP4zne+gxMnTmBhYYHrAWnT/Pw85ubm0NbWBgDo6OjAgwcPNvfLnQQCgc31pi6XC319fXj8+PHmC6poNIrl5WVcunQJFy5c4NmezvHTGahqFosFPp8PmUwGTqcTkiRheHgYExMTUFUVqVQK3/72t3Hy5Mm63aGls7MTnZ2dm+sBP/74Y64HNLGtZ3sbRFHcPOv7wQ9+sOPjfD4ftlYfv9+PeDyOp0+fwul0wuv14vd+7/f2dKmetINnfLQnXV1d24L/xqcwPHjwAL//+7+P06dPN+S2ZK2trThx4gQuXryIlpYWfPbZZ/j000+xuLgIZmzzeP5sb0Opsz6Px/PCftLT04OpqSm0t7fjxz/+MYeegXDw0Z5svUS0cWnz+PHjePPNN7GystLw7WEHNK+dzvY2bD3r20lLSwskSdp2aTMWi+HP/uzP4Ha7kcvV5xMhqDl4qZP2ZOOVdTKZxMrKyualTUVRcO3aNczOzmLfvn0N366NDhgMBjfXA05MTHA9oIFtnO1t3CHoeRtnfadPn37hHZ6iKMLv9yOVSiGRSGy7tDkzM4MvvvgC58+f540UDIJnfLQnHo8HhUIBmUwGP/7xjzcvbVosFrz88su4e/cukslkU7eR6wGNb7ezvQ2lzvoCgQBCoRBGRka2Xdrcv38/vF4v7t69W5dtp8bjOj7aE1VVcefOHRw4cAAtLS0v/P7Gurtz585p5tUy1wMaz9zcHN5///1tZ3uyLOPmzZt4+eWXN39NUZTNe8U+f9YXiUQQj8cxNDT0wr5QKBTwq1/9CsPDwwgGg/X9w1DdcfBR3d24cQMWiwWjo6PN3pRteF9Q43j//fcxNTWF1tbWzV+TZRn379/HiRMntn3t6uoqTp8+XfQdnsXE43GMj49jbGxs2/ch/eHgo7orFAq4du0aRkZGmtL7ysH7gurbxMTEC2s4FUXB9evXce7cuRe+vqWlBQcPHqz4+8zMzGBqaoq9T+c4+Kgh9PJqmZ8PaByyLOPq1atFPwmkWjdu3IAkSTh58mRNn5cah29uoYbweDw4fPgwvvzyS00vK+B6QCpldHQUKysrCIVCzd4UqhIHHzVMf38/3G437t271+xNKYnrAakYi8WCM2fO4N69e01/xzJVh4OPGmp0dBRLS0uYnZ1t9qaUhfcFpZ14PB4cOXIEX3zxBV8I6RAHHzWUltb3VYrrAWmr/fv3w+PxcH2fDnHwUcPppfcVww5IG9j79ImDj5pCT72vGHZAYu/TJw4+ahq99b5i2AHNjb1Pfzj4qGn03PuKYQc0J/Y+feHgo6bSe+8rhh3QfNj79IODj5rOCL2vGHZA82Dv0w8OPtIEo/S+YtgBzYG9Tx84+EgTjNj7imEHNDb2Pu3j4CPNMGrvK4Yd0LjY+7SNg480xci9rxh2QONh79M2Dj7SHKP3vmLYAY2FvU+7OPhIc8zU+4phBzQG9j5t4uAjTTJb7yuGHVD/2Pu0h4OPNMuMva8YdkD9Yu/THg4+0jSz9r5i2AH1ib1PWzj4SNPY+4pjB9QX9j7t4OAjzWPv2x07oH6w92kDBx/pAntfaeyA2sfepw0cfKQb7H3lYQfUNva+5uPgI91g76scO6A2sfc1Fwcf6Qp7X3XYAbWHva95OPhId9j7qscOqB3sfc3DwUe6xN63N+yA2sDe1xwcfKRL7H21ww7YXOx9jcfBR7rF3ldb7IDNw97XWBx8pGvsfbXHDth47H2NxcFHusfeVx/sgI3F3tc4HHyke+x99ccO2BjsfY3BwUeGwN7XGOyA9cfeV38cfGQY7H2Nww5YP+x99cfBR4bC3tdY7ID1wd5XXxx8ZCjsfc3DDlhb7H31w8FHhsPe11zsgLXD3lcfHHxkSOx9zccOuHfsffXBwUeGxd6nDeyAe8PeV3scfGRY7H3aww5YHfa+2uLgI0Nj79MmdsDKsffVDgcfGR57n3axA5aPva92OPjIFNj7tI0dsDzsfbXBwUemwN6nH+yAu2Pv2zsOPjIN9j59YQcsjr1vbzj4yFTY+/SHHfBF7H17w8FHpsPep0/sgNux91WPg49Mh71P/9gB17H3VYeDj0yJvc8Y2AHZ+6ohqGbZO4h2cOPGDVgsFoyOjjZ7U6gGFEXB3NwcJicnkc/nEQqF8Pbbb0OSpGZvWl3F43GMj49jbGwMra2tzd4czeMZH5kae5+xPN8Bo9GoKToge19lOPjI1Nj7jKuzsxMjIyOm6YDsfeXj4CPTY+8zNjN1QPa+8nDwEYHr+8zADOsBub6vPBx8RL/D3mcORl8PyN5XGgcf0e+w95mPUdcDsvftjoOPaAv2PnMyYgdk7yuOg4/oOex95mWkDsjeVxwHH9EO2PvMzSgdkL1vZxx8RDtg76MNeu+A7H0v4uAjKoK9j7bScwdk79uOg49oF+x99Dw9dkD2vu04+IhKYO+jnRTrgBMTE5rsgOx93+DgIyqBvY9K2doB8/m8Zjvg/v374fV6Td/7OPiIysDeR+XYqQOOj49rqgOeOHHC9L2Pg4+oTOx9VK6tHbCvr09THZC9j4OPqCLsfVQJrXZAs/c+Dj6iCrD3UbW01gHN3Ps4+IgqxN5He6GlDmjW3sfBR1SFjd5nxlfLVBta6IBm7X0cfERVGh0dxfLyMnsf7UmzO6AZex8HH1GV2Puo1prVAc3W+zj4iPaAvY/qoRkd0Ey9j4OPaI/Y+6heGtkBzdT7OPiIaoC9j+qpUR3QLL1PV4NPUVRk8jLSORmZvAxF0cYtgIjY+6hR6t0BK+l9ej0mC6pWbiBXRCyVQyiaRng1jUgyB3nL5kqCAH+rDQGvE0GfE20uWxO3lAiYnp7G06dPce7cOUiS1OzNMTVZlnH16lVcvny52ZtSV7lcbnO/a21txYEDB+D3+yEIQtXPWSgU8Ktf/QrDw8MIBoPbfs8Ix2TNDr751TTuhFaxEM9AEgW47BY4rRIk8Zu/TFlRkc7LSGULkBUV3R4HTgS96PE6m7jlZHY3btyAJEk4efJkszfF1Mwy+DYoioK5uTlMTk5CURQMDQ0hGAxW/QIsHo9jfHwcY2NjaG1tNdQxWXODL5OXcfNZFI/CSbidFrTaLWW9clFVFclsAfF0HocCbpzu88Fh5StuarxCoYBr165hZGQE+/bta/bmmJbZBt9Wy8vLmJycRDQaRX9/PwYHB2G32yt+nmfPnmEpGkOmtcdQx2RNDb5YKoePJhaRzSvocNsgVnGqrqgqlhM52K0iLh7u0uypNhnb86+WqfHMPPg2rK2tYWpqCrOzswgEAhgaGoLH4yn78UY9JmvmzS2xVA5X7y9AEAC/x17Vf2AAEAUBfo8dggBcvb+AWCpX4y0lKo3r+0gLWlpaql4PaORjsiYGXyYv46OJRdgkAW6HtSbP6XZYYZMEfDSxiEyeBx5qPK7vI62odD2g0Y/Jmhh8N59Fkc0rNfsPvMHtsCKbV3DzWbSmz0tULq7vIy3Zuh5wdHS06HpAox+Tmz745lfTeBROosNdn+u+HW4bHoYTCK9mKnpcLpfDvXv3sLKyUpftInPg+j7Sqo6ODrz66qs4d+7ctvWAM0sJTR6Ta6npg+9OaBVup6Xq68eliIIAj9OK26FYWV+fy+Vw584d/MM//AN+/vOfIxrl2SLtDXsfadnWDuj1enFnNqapY3I9WJr2nbEeTxfiGQS8jl2/7sb1a8hm118diKIIq9WGVq8XPcF+uL1tJb9Pq92C8GoGsVSu6DuKcrkcHjx4gM8//xy5XA5+vx/5fL7iPxPRTvr7+7G8vIy7d+9yfR9pks1mg697H1bCcxUdkwUBsFht8Pra0X9wBDZb6WUT5RyT66mpgy8UTUMShbLvMODr6ITFakNiNYblxQWsRBZw8OgoOrsCuz5OEARIooBQNP3Cf+SdBl41612IShkdHcW1a9cwOzvL9X2kSdUck+0OJ5YjC1haCAMAho+OlnzcbsfkRmjq4AuvpuGyl78JXT1BtPu7oKoqHn99B0sLYTx5cB++Dn/JuxO47BaEV9M4vs8LgAOPGm+j942Pj8Pr9XJ9H2lOtcfkFrcHkxP3sJYo/36hzx+TG6lpg09RVESSObS3VD7tBUFAsP8AlhbCKBQKSKxG0dbeuetjnFYJkWQOmUwWjx495MCjptja+3g/T9KSao/JiqJgLbk+8Fyt7rIft3FMVhQVolifnlhM0wZfTlYgq+q2+7xVwub45hp0oYwWJ4kCnoVm8dfv/x/I+Sw6OjpgtVoxPT1d9DGLi4v46quvsLy8XNU2EhUzPT2NZ8+e4dChQ83eFMOSZRkPHjxAZ+fuL4ppnSJIkNXOio7JD+7e2vyxp82HweEjZT9WEgXIqoqcrMAhNvYFYNMG315vlJbLfPNWWIu1vLUmPp8Pw2dfw6OJexBFER0dHbu+4pZlGfv27UN/f//eNpboOfv27cOnn34Km82Gnp6eZm+OIcmyjKmpKf77LVNeFfB4Ol3RY3wdnVBUFasry0jGY0in1mC1VXbG2IybZjZt8O3lnbKqqiI0PQlgvZu4vb6yHud0OfGjc5fw+ve/izt37uCrr76CIAjo6uracQCura2hs7MT3d3d1W8sURE//OEPMT4+joMHD7L31YEsy/D5fPz3W6ZMXgamQxU9ZqPxPf76LiLhOTx9PIHRl1+r6DnqtGpiV01bx2eTREiCALmCDy5cnA/h8dd3ceuz32BpIQxBAAZHjpbVSWRFhSQIsEkiWlpacPbsWfzJn/zJ5t0L5ufnucaKGorr+0hLqjkmb+gbPABBANYSCcRWlsp6zNZjcqM1bfCJ4voHFqYruGdbdHkJy4thqIqCjq5uHH/pWyWXMmxI52X4W23bIioHIDUb7+dJWlHNMXmD3eFEZ3cvAGB2+klZj9npmNwoTV3OEPA6cTsUQ2uJt8++9O0Le/5eqWwBB/07X07aGIAnTpzAnTt3cOvWLaiqimw2W/IO5kR7xfV9pBV7OSYfPHIcB48cL/t77XZMrremDr6gz4mbM1Goqlr2gslqqKqKJ09n4M+EkJm3I5fLIZvNIpfLIZfLIZ/Pb/txJpPB3NwcZmZmcPr0aRw8eLBu20bE9X3UbIVCATMzM3j6eBop+z6obnvdj8myoiLoa84nszd18LW5bOj2OJDMFmp+F/CtktkChGwc//x/30EgEEBPTw+8Xi8kSYIoipAkCZIkobW1FaIowu/3o7+/H319fXjppZfqtl1EG7i+j5ohnU7jyZMnmJmZgd/vx9grp/DbuWxDjsndHkfTPpS26TepPhH0IpEuQKnTJUVFVRFP5/Gf3/wefvrTn8Lj8WBtbQ2zs7PI5XJoaWmBy+WCw+GA1WqFJEkQBAFLS0s4c+YM3O7yF2QS7QV7HzVKLBbDl19+if/4j/+Aqqq4cOECzpw5A5/P17Bj8miwrS7PX46mnvEBQI/XieFAK55GUvB7an8HleVEDocCbvT6XOj43vcQj8eRSqVgtVoxPz+PUCgEv9+P7u5uWH+3HlBRFBQKBRw7dqzm20O0G/Y+qhdVVREOhzE1NYVMJoPBwUGcPHkSFsv2MdCoY3KpG2HXU9MHHwCc7vNhPpZBIpOv6el1IpOH3SridN/6Oj+73Y5Lly7hnXfegSiKGBkZQTqdRjgcxu3bt9He3o5AIIBUKoWhoSF4vY2/hxyZG3sf1dpGv3vy5AnsdjsOHDiAQCCwa8Nr1DG5WZp+qRMAHFYJFw93ISerSGRq81FAiUweOVnFxcNdcFi/6SUejweXL1/GysoKcrkcnE4nBgcHMTo6CrvdjomJCdy9e5evtqlpuL6PaiGdTuP+/fv45S9/iWg0ipdeegnnzp1DT09PyTeuNPKY3AyaGHzA+htd3jjaDVUFIvFs1deXFVVFJJ6FqgJvHO3eMZ4GAgFcvHgR8/PzUBQFAGC1WtHb24uhoSEcPHgQS0tL+OSTTzAzM7P5NUSNwt5H1dqt31XCZQG6cvMNOSY3miYudW5oc9lw+UQPbj6L4mE4AY/Tila7pay31aqqimS2gHg6j0MBN073+XZ9VXHkyBGsrKzgq6++Ql9f3+avx+Nx/OhHP8Lg4CCWlpYwOTmJiYkJ9Pf3Y2BggJ/kQA3D3kflKrfflSOXy+Hhw4e4cuUKEokE/utf/zfcDMXqfkxuJE0NPmD9FPu1oU4MdrTidiiG8GoGkijAZbfAaZW23TlcVlSk8zJS2QJkRUW3x4HXhjrLjqZnz55FLBbD/Pw8uru7kclk0NLSsjkIOzs70dnZiWQyiampKXz88cfo6enB0NAQ3+1JdcfeR6VU0++K2Rh4n3/+OaLRKJaXl/Hqq6/CYbM07JjcKIKq8VuTxFI5hKJphFfTiCRzkLdsriSs32In4HUi6HNWdQqdTqfx3nvvQZZlpFIpjI2N4cSJEzt+bS6Xw/T0NJ48eQKPx4MDBw7A7/dX/WcjKsf09DSePn3K9X0VkmUZV69exeXLl5u9KTX3/Pq7oaGhii9lbtg68LLZLHw+Hx4/fgyn04nz58/j1Vdf3fb19T4mN4LmB99WirL+2U2qun5Hb5sk1uQ+b9FoFD/72c+gqir+9E//FA7H7q9OFEXB7OwspqamoKoqhoaGEAwGIYqaSaZkMDdu3IAkSTh58mSzN0U3jDj4YrEYJicnEYlE0NfXh8HBQbhcrqqe6/mBt/Gh3I8ePYLdbocoirh06RIGBweLPke9jsn1prlLnbsRRaEuH1jo8/k23+lZauitb4eIvr4+9PX1sQNSQ7D3mVct+x1QfOABQDgcRi6Xw4EDBzA/Pw+Px7Prc9XrmFxvuhp89RQMBhEMBit+HDsgNQJ7n/nUst9t+PrrrzE+Pv7CwAOAZDKJ+fl5HD16dPN7lBp8esXBVyOtra0YHR3F4cOHMT09jfHxcXZAqinez9Mcnu93L730UtX97nnZbBapVAqdnZ3bhl6hUMDk5OTmFatMJgOPx7N5Nyuj4eCrMZvNhuHhYRw4cACzs7O4f/8+OyDVTH9/P5aXl3H37l32PoN5vt9duHCh6n5XzKlTp9DZ2YkPP/wQa2tr6OrqAgA8efIEPp9vc8Cm0+mqroDpBQdfnbADUr2w9xlHrftdOYLBIH7yk5/go48+wszMDCRJ2ux6G1KpFAKB8j7kW484+BqAHZBqib1P/+rR7yrR0tKCy5cv49q1a3j33XfxyiuvvHA1qq2trSHb0gy87tZAGx3w4sWLcLlcGB8fx6effopIJNLsTSOd4f089Wkv98+sNUVRoCgK/vzP/xyKomBxcXHb7xv1jS0Az/iagh2QaoG9Tz8a0e8qdevWLQQCARw7dgwjIyOblz4DgQAsFgtaWlqaun31xMHXROyAtFfsfdrVjH5Xro1tOnPmDIBvLn3eunULv/71rzE4ONjwM9BGav7fAAFgB6TqsPdpT7P7XSmxWAyPHz/GuXPntl1dkiQJZ86cQSAQQD5fm48i0ipd3bLMTHhfUKoE7+f5okbfsqyW98+sl3w+j2vXruHYsWOGftdmKTzj0yh2QKoEe1/zaLHfFbPR9cw89AAOPs1jB6Rysfc1jpb7XTHPdz0z0+7fEr2AHZB2w95Xf1rvd8UU63pmxcanY+yAtBP2vnW1bHx66HfFsOu9iGd8OsYOSDth76sdPfW7Ytj1XsTBZwDsgPQ89r7q6bHfFcOutzP9/U3SrtgBCWDvq4Ze+10x7HrFsfEZHDuguZm595Xb+PTc74ph19sdz/gMjh3Q3Nj7ijNCvyuGXW93HHwmwQ5oXux93zBSvyuGXa804/xtU9nYAc2Fvc94/a4Ydr3ysPERO6BJmK33ybKMDz74ACMjI4bqd8Ww65WPg482KYqC2dlZTE1NsQMa1I0bN2CxWDA6OtrsTamrWCyGhw8f4t///d/xh3/4hxgcHDRMvyvmt7/9LVwuF44dO9bsTdE8Dj7a0UYHXF1dZQc0kEKhgGvXrmFkZMRwve/5frd//35MTEzgzTffbPam1d3U1BRmZ2cxNjbGF6plYOOjHbEDGpMRe1+xfqcoCh49etTszas7dr3K8YyPysIOaCxG6H2l1t81+vP4moFdrzo846OycD2gsWys77t3757uep+R199Viuv1qsPBRxXhekDj0NP6PjOsv6sU1+tVz7x7De0ZO6C+6aH3mWX9XaXY9faGjY9qhh1Qn7TY+2px/0yjNj52vb3jGR/VDDugPmmp97Hflcaut3ccfFRz7ID608zex35XPna92uCeRXXFDqgPzeh97HeVYderHTY+aih2QG1rRO9r1OffGanxsevVFs/4qKHYAbWtnr2P/a567Hq1xcFHTcEOqF217H3sd3vHrld73Puo6dgBtaUWvY/9rjbY9eqDjY80hx1QG6rpfY3qd+XQe+Nj16sfnvGR5rADakMlvY/9rvbY9eqHg480ix2w+Xbrfex39cOuV1/cQ0kX2AGbY6fex35XX+x69cfGR7rEDthY09PTePDgAXp6ejA7O9v0flcOPTY+dr3G4Bkf6RI7YOPEYjEsLS3h4cOHSCaTuHTpEvtdnbDrNQYHH+kaO2B97NTv/uIv/gLXr19HNBrl4KsDdr3G4eAjw2AH3LtS/U7rn9+nV7FYDI8ePcL58+d5taIB2PjIsNgBy1fJ+jstfn7fTvTS+Nj1Go9nfGRY7IClVbP+Tkuf32cE7HqNx8FHhscOuF0t1t818/P7jIRdrzk4+MhUzNwBa7n+rhmf32c07HrNw8ZHpmaGDljP+2dqufdpufGx6zUXz/jI1IzcARtx/0z2vuqw6zUXBx8RjNMBm3H/TPa+yrDrNR8HH9Fz9NgBm3n/TPa+8rHraQMbH1EJWu6AWvr8O631Pq01PnY97eAZH1EJWuyAWvz8O/a+3bHraQcHH1GZmt0B9fD5d+x9O2PX0xbt/Ish0pFGdkA9ff4de9+L2PW0h42PqAbq0QG11O8qpYXep4XGx66nTTzjI6qBWnZALfa7SrH3rWPX0yYOPqIaKtYBBwYG0N/fX7QD6qHfVcrsvY9dT7v0+6+KSOPK6YB66neVMnPvY9fTNjY+ogbZ2gHtdjskSUIymdRdv6tUs3pfsxofu5728aUIUYPYbDb4/X74fD7MzMzg8ePHAAC/3w+v19vkrauf/v5+uN1u3Lt3r9mb0hDsetrHS51EdbZTv/ujP/ojWCyWijqgnpml97Hr6QMHH1GdlNPv9Hhf0GqYofex6+kHGx9Rje1l/Z2W7wtaC43sfY1sfOx6+sIzPqIaqcX6Oy3eF7SWjLq+j11PXzj4iPagXuvvql0PqAdG633sevrDwUdUha39zuFwYGhoqG7r74zWAY3U+9j19ImNj6gCWrh/plE6YL17X70bH7uefvGMj6gMWrp/plE6oN57H7uefnHwERWh9ftnGqED6rX3sevpmzb+BRNpSCP7Xa3otQPqsfex6+kfGx/R72ih39WK3jpgPXpfPRofu54x8IyPTE9L/a5W9NYB9dL72PWMgYOPTEnr/a5W9NQBtd772PWMw1j/yolK0GO/qxWtd0At9z52PWNh4yNTMFK/qxWtdsBa9b5aNT52PePhGR8ZmhH7Xa1otQNqrfex6xkPBx8Zjln6Xa1osQNqpfex6xkTjwRkGGbud7WilQ6ohd7HrmdcbHyke+x39dPsDriX3reXxseuZ2w84yPdYr+rv2Z3wGb1PnY9Y+PgI11hv2uOZnbARvc+dj3j49GCGiqTycDhcFT8OPY77Wh0B2xk72PXMwc2PmqYx48f4/r16/jJT35S9vBjv9O+RnXASntfpY2PXc88eMZHDRGJRPDhhx8in89jcnISx44d2/Xro9Eopqam2O90oFEdsN69j13PPDj4aBtFUZGTFagqIAiATRIhinu7nLi2toZf/OIX8Hq9sFqt+OKLL3D48OEXXrWz3+lbIzpgub1PUVTkCgpkQUImL5fcj9n1zIWXOgmxVA6haBrh1TQiyRzkLbuEJAjwt9oQ8DoR9DnR5rJV9Nz5fB5XrlxBLBbbvPz17NkzXL58GQMDAwDY74xsowPOzc3VrAPG43GMj49jbGxsW+/baT+WZRmSJO26H8diMXz22Wc4f/48ryqYBAefic2vpnEntIqFeAaSKMBlt8BplSBteWUsKyrSeRmpbAGyoqLb48CJoBc9XmfJ51dVFZ988gkePHiw7dV5PB6H3W7H5cuX8fTpU/Y7E6h1B9za+xaTuar3406XhV3PhDj4TCiTl3HzWRSPwkm4nRa02i1lnV2pqopktoB4Oo9DATdO9/ngsBZ/k8GtW7fw61//Gvv379/2/MlkErdu3cKpU6cwOjqKwcFBvtI2CUVRMDs7i6mpqT13wM+/vImI4sJy3l71fuxRUxjyAieP796cyVg4+Ewmlsrho4lFZPMKOtw2iFVcTlRUFcuJHOxWERcPd+14+XN6ehpXrlzBvn37YLFYoKoqotEowuEw8vk8bDYbRkdHcenSpVr8sUiHNjrg6upqxR2wUfsxGRMHn4nEUjlcvb8AmyTA7bDu+fkSmTxysoo3jnZvO2isrKzgnXfeQVtbG2w2GyKRCBYWFmC1WhEIBODz+aCqKubm5vDHf/zH8Hg8e94W0q9KO2Cj9mMyLg4+k8jkZfzizjwEATU5WGxIZPJQVeDyiR44rBLS6TTee+89ZDIZZDIZRCIReL1eBAKBFxYez8/PY3R0FGfPnq3Z9pB+ldMBG7Ufk7Fx8JnE+NQSnkZS8Htqf2upSDyLAb8Lr/b78C//8i+4f/8+JEmC3+9Hd3d30ctXhUIBS0tLePvtt6u6mwsZ024dsBH78WtDnTV/btIWDj4TmF9N48N7Cwi0OapqIaUoqor5WBot0Ulcu/oBurq60N7evrlOb+supqoqJEmCKIoQBAFTU1P4zne+g9dff73m20X6t7UDurv78GjNWff9+I1jPQh4+ULMyLgy2ATuhFbhdlrqcrAAAFEQ4HFa4fQcxd/8zVkoioJCoQBZll/4/3w+j0wmg5mZGdy6dQvRaBShUKgu20X6t/W+oP/v88dw+9x1349vh2IIeLm0wcg4+AwulsphIZ4p6xXsjevXkM1mXvj10VdeQ0vr7ouOW+0WhFczKIg2tLXu/AYBRVE2B97KygqGh4cxMDCA3t7e8v4wZFoF0QY429BqL/+QtXV/PvXqGJwtLSUfs7Efx1I5vtHFwDj4DC4UTUMShYruguLr6ITD+c26Oqu19JsIBEGAJAoIRdMvHDA2Bt6nn36KlZUV+Hw+9PX1AVi/JydRKZXux/HYyrYXcZGFOewfGi75uN32YzIODj6DC6+m4argVTIAdPUE0e7vqvh7uewWhFfTOL7PC2D3gUdUiUr340h4HgDQ4nZjLZHAUni+rMEHvLgfk/Fw8BmYoqiIJHNob6nslevifAjx2MrmzweGD5f1OKdVQiSZQ6EgIxR6xoFHNVHpfqwoClYiCwCA/gOH8PDuV8hmM4jHVuBpay/5+I39WFHUPd+gnbSJg8/AcrICWVW33bOwHNHlpW0/L3fwSaKARDKJ//4//ieiS4vw+XxwuVyYn58v+phEIoFEIgGvl6+uaWc5WUU4IsB/cH9ZXx9djqBQKMBqs8HT1g5fhx+RhXlEwvNlDT5JFCCr659S4hC5ps+IOPgMrNqFKiPHT1V1qRMAHA4nvnX2LO7euoFsNgu3273rbaisViv8fj/PCKmoTF6BO5Mo++uXwnMAAF+HH4IgoN3fjcjCPFYiCxg8dKTs+4JyoZdxcfAZWDM+1UeURJw/dw7fPfdtPHz4EJ999hnW1tbg9/uLDsD29nbeGZ+KyuRlOObyZX1tIZ9HbGX9isXi/CwW52e/+b1CAdHlCDr83WU9Fz8Vy7g4+AzMJomQBAGyUtnlzucbn79nX8nlDMD6R79IgvC7D/2UcOzYMRw6dGhzAGaz2V0HINFOKtmPlxfDUJT1myR4fd9c1kytrSGTTmEpPFdy8G3dj8mYOPgMTBTXP3wzmZMrWv/0fOPztLWXNfjSeRn+Vtu2NwRYrVYOQNqTSvbjyMJ6T+7uDaL/4Mjmr8djK7h38wvEVpaQz+VgtRV/o8xO+zEZC29ZZnB3Z1dxOxRDl6f+t2BajGcwGmzb9W3g+Xx+2wC0WCwYHBzED37wg7pvH+mX1vZj0jeeyxtc0OeErKio9+sbVVUhKyqCvt0/mX3jDPDtt9/GhQsXIEkSLBZeeKDdaW0/Jn3jEcfg2lw2dHscSGYLNf0Yl+clswV0exxl3+1i6yXQfL68Ny6QeWl1PyZ94hmfCZwIepFIF6DU6dWyoqqIp/MYDbZV/Fir1QqXy1X6C8n0tLwfk75w8JlAj9eJ4UArlhO5ujz/ciKHQwE3P8qF6or7MdUKB59JnO7zwW4VkcjU9rJiIpOH3SridJ+vps9LtBPux1QLHHwm4bBKuHi4CzlZrdlBI5HJIyeruHi4Cw4rb+1E9cf9mGqBg89E2lw2vHG0G6oKROLZqluJoqqIxLNQVeCNo918IwA1FPdj2iuu4zOhTF7GzWdRPAwn4HFa0Wq3lPU5Z6qqIpktIJ7O41DAjdN9Pr5CpqbhfkzV4uAzsfBqBrdDMSzEM5BEAS67BU6rtO22ULKiIp2XkcoWICsquj0OjAbb+AYA0gzux1QpDj5CLJVDKJpGeDWNSDIHecsuIQnrt4sKeJ0I+py8HESaxf2YysXBR9soyvrnkKnq+t3p1284zXsWkr5wP6bdcPAREZGp8F2dRERkKhx8RERkKhx8RERkKhx8RERkKhx8RERkKhx8RERkKhx8RERkKhx8RERkKhx8RERkKhx8RERkKhx8RERkKhx8RERkKhx8RERkKhx8RERkKhx8RERkKv8fXekMir3ol5kAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from pgmpy.estimators import TreeSearch\n",
    "\n",
    "# learn graph structure\n",
    "est = TreeSearch(df_data, root_node=\"R\") # Try n_jobs=1 if estimate() fails below.\n",
    "dag = est.estimate(estimator_type=\"tan\", class_node=\"A\")\n",
    "nx.draw_circular(\n",
    "    dag, with_labels=True, arrowsize=30, node_size=800, alpha=0.3, font_weight=\"bold\"\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## To parameterize the learned graph from data, check out the other tutorials for more info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<TabularCPD representing P(R:3 | A:2) at 0x7f43c56f1820>,\n",
       " <TabularCPD representing P(C:2 | A:2, R:3) at 0x7f43c57328e0>,\n",
       " <TabularCPD representing P(B:3 | A:2, R:3) at 0x7f43c56be4c0>,\n",
       " <TabularCPD representing P(D:3 | A:2, R:3) at 0x7f43c56f1b20>,\n",
       " <TabularCPD representing P(E:2 | A:2, R:3) at 0x7f43c56becd0>,\n",
       " <TabularCPD representing P(A:2) at 0x7f43c56f1a00>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from pgmpy.estimators import BayesianEstimator\n",
    "\n",
    "# there are many choices of parametrization, here is one example\n",
    "model = DiscreteBayesianNetwork(dag.edges())\n",
    "model.fit(\n",
    "    df_data, estimator=BayesianEstimator, prior_type=\"dirichlet\", pseudo_counts=0.1\n",
    ")\n",
    "model.get_cpds()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
